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Abstract 

A generalized  Monte  Carlo  Analysis  Program  (i-IGAP)  has  been 
developed  for  linear  or  extended  Kalman  filter  design.  The  computer 
program  is  similar  in  structure  to  the  Air  Force  Avionics  Laboratory 
developed  General  Covariance  Analysis  Program  (GCAP)  so  users  can 
easily  transition  from  a covariance  analysis  to  a Monte  Carlo  analysis 
of  a (extended)  Kalman  filter,  A detailed  users*  guide  for  MCAP  is 
Included  in  Appendix  A.  In  addition,  this  study  treats  the  high  accur- 
acy tracking  of  a satellite  from  an  aircraft.  The  purpose  of  the  study 
is  to  evaluate  the  feasibility  of  a reduced  order  system  model  for 
implementation  in  an  extended  Kalman  filter  formulation  whose  estimates 
vfould  be  used  to  aid  the  tracker.  The  six  state  model  accomplishes 
tracker  state  estimation  by  exploiting  the  information  already  avail- 
able in  the  tracking  geometry,  dominant  modes  of  satellite  dynamics, 
and  the  range  measurement.  Tracker  state  estimation  is  accomplished 
in  the  line-of-sight  coordinate  frame.  A Monte  Carlo  analysis  was 
performed,  evaluating  the  filter  against  a 42-state  truth  model.  With 
some  proposed  jriOdificatu.ons,  it  was  concluded  what  the  six  sv^at^,  filter 
is  feasible  and  warrants  further  study. 
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A GEI^EIULIZSD  tlC*NTE  CARLO  AI'JALYSIS  PROGRAl'I  FOR  KALIWI  FILTER  DESIGN 
WITH  APPLICATION  TO  AN  AIRCRAFT-TO -SATELLITE  TitAGKIKG  FILTER 

I.  Introduction 


Overvi ew 

A large. class  of  estimation  problems  is  concerned  with  finding  an 
optimal  estimate  of  some  quantity  (an  unknown  parameter,  a random 
variable,  or  a random  signal)  given  noise-corrupted  measurements  of 
a funcd^ion  of  this  quantity  (corrupted  by  noise).  In  particular,  there 

Znany  problems  with  aerospace  applications,  such  as  navigation, 

Lance,  and  weapon  delivery  systems,  which  require  that  noise- 
corrupted  measur-ements  be  used  to  estimate  certain  physical  parameters 
precisely.  For  instance,  the  techniques  used  to  combine  these  external 
measurements  with  Inertial  Navigation  System  (iNS)  outputs  fall  into 
two  general  categories:  conventional  continuous-feedback  damping  and 
Kalman  filter  damping  (Ref  l:l).  The  trend  in  recent  years  has  been 
toward  extensive  use  of  Kalman  filter  techniques. 

In  complex  aerospace  systems,  the  number  of  parameters  needed  to 
describe  the  system  (the  true  system  model)  accurately  may  be  extremely 
large.  The  real  time  implementation  of  this  true  system  model  within 
the  "optimal"  Kalman  filter  is  often  not  practical  becai:ise  of  the 
large  number  of  parameters,  memory  requirements,  and  resultant  compu- 
tational b\jrden  imposed  on  the  airborne  computer.  In  oirder  to  obtain  a 
Kalman  filter  which  is  computationally  feasible.  Intentional  modeling 
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approximations  are  introduced  by  deleting  nondominant  states  and 
terms  from  the  system  model.  Thus  a "suboptimal"  Kalman  filter  based 
on  this  reduced  order  model  will  give  poorer  performance,  from  a 
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statistical  standpoint,  than  the  optimal  filter.  Computer  simula- 
tion techniques  are  used  to  develop  the  required  compensation  (filter 
tuning)  and  to  study  the  effects  of  uncertainties  (an  error  analysis) 

In  the  true  system  model.  Two  types  of  siraulation  techniques  are 
commonly  used!  the  covariance  analysis  and  the  ilonte  Carlo  analysis. 

A covariance  analysis  generates  the  second  order  statistics  of 
the  error  between  suboptimal  Kalman  filter  state  estimates  and  the 
corresponding  parameters  from  the  true  system  model.  A General 
Covariance  Analysis  Program  (GCAP)  has  been  developed  by  the  United 
States  Air  Force  Avionics  Laboratory  and  is  coming  into  widespread 
use  (P.ef  2).  The  covariance  analysis  is  somewh-at  limited  because 
stringent  assumptions  (see  Appendix  A)  are  necessary  for  the  analysis 
results  to  be  a valid  depiction  of  the  error  characteristics. 

The  Monte  Carlo  analysis,  on  the  other  hand,  actually  conducts  a 

sample-by-sample  simulation  using  random  number  generators  and 

shaping  filters  to  generate  the  random  error  sources.  The  nqiay  system 

measurements  are  then  processed  by  the  suboptimal  (extei^ded)  Kalman^ 

! 

filter  algorithm  to  generate  the  filter  estimates.  If  mWiyJof-'ihes e 
simulations  are  conducted,  then  the  statistics  of  the  error  between 
the  filter  estimate  and  the  truth  model  can  be  computed.  Because  of 
the  larger  number  of  simulations  necessary  to  determine  the  performance 
of  cne  set  of  Kalman  filter  parameters,  a lionte  Carlo  analysis  can 
require  the  expense  of  large  amoimts  of  computer  resources. 
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Consequently  when  a covariance  analysis  can  be  performed,  the  Monte 
Carlo  analysis  is  viewed  as  a step  to  be  performed  after  a covariance 
analysis.  The  Monte  Carlo  analysis  is  thus  used  to  fine  tune  the 
Kalman  filter  and  validate  the  filter  performance  after  the  covariance 
analysis  has  been  performed. 

The  extended  Kalman  filter  cannot  be  evaluated  or  tuned  totally 
on  the  basis  of  a covariance  analysis  (like  GCAP)  because  the  gain  and 
error  covariance  matrices  depend  on  the  time  history  of  the  state 
estimates  (which  are  not  available  in  a covariance  analysis).  Never- 
theless, the  covariance  analysis  is  used  to  tune  a linearized  Kalman 
filter  operating  over  a nc  '.i  nal  state  trajectory  as  an  approximation 
to  the  extended  Kalman  filter  with  small  deviations  from  this  nominal. 
Subsequently,  the  Monte  Carlo  analysis  is  conducted  to  investigate 
filter  performance  thoroughly  (Ref  3‘9-3S). 

Since  a Monte  Carlo  analysis  of  promising  filter  designs  should  be 
performed  after  the  covariance  analysis,  it  is  desirable  to  have  avail- 
able a general  Monte  Carlo  Analysis  Program  (lIGAP).  This  program 
should  be  similar  in  structure  to  GCAP  so  that  users  can  easily 
transition  from  a covariance  analysis  to  a Monte  Carlo  analysis  of  the 
Kalman  filter  under  study.  Also,  if  possible,  problem  dependent 
computer  subroutines  should  be  applicable  to  both  GCAP.  and  MCAP.  The 
development  of  this  general  Monte  Carlo  Analysis  Program  is  the  primary 
objective  of  this  study. 


Ob.iectlvcs  of  the  Study 

The  primary  objective  of  this  study  is  to  develop  a general  Monte 
Carlo  analysis  program  -which  can  be  used  to  evaluate  llalman  filter 
and  extended  Kalman  filter  performance.  In  order  -to  demonstrate 
the  validity  of  the  simulationi  the  following  second2Lry  objectives  will 
be  performed  to  investigate  a specific  problem. 

1.  To  develop  a simulation  of  the  aircraft  and  satellite  dynamics 
using  the  truth  model  developed  by  Mitchell  and  modified  by 
Mann. 

2.  To  perform  a Monte  Carlo  performance  analysis  of  the  six 
state  reduced  order  filter  investigated  by  Mann. 

3.  To  analyze  the  performance  of  the  reduced  order  filter. 

Problem  Statement 

The  high  accvaracy  tracking  of  one  accelerating  vehicle  from 
another  vehicle  has  many  military  applications.  These  include*  aircraft 
to  aircraft  tracking,  aircraft  to  missile  tracking,  aircraft  to  satel- 
lite tracking,  missile  to  aircraft  tracking,  etc.  This  project  will 
consider  the  problem  of  tracking  a satellite  from  an  aircraft.  Tracking 
is  usually  accomplished  by  measuring  range  and  angle  information  with 
a steerable  radar.  These  measurements  eire  then  used  to  estimate  the 
state  (position,  velocity,  acceleration,  or  equivalent  parameters) 
from  the  aircraft. 

Since  the  radar  observations  of  the  satellite  are  uncertain  due 
to  noisy  measurement  da-ta,  the  state  estimates  ere  not  exact  (determin- 
istic) quantities.  If  the  measurement  corruptions  are  described 
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statistically  and  the  system  dynamics  axe  described  by  a linear  mathe- 
matical model,  then  a Kalman  filter  will  provide  the  si  wistics  of  the 
state  estimates.  However,  since  the  problem  to  be  investigated  is  non- 
linear, the  extended  Kalman  filter  formulation  has  been  chosen  to  deter- 
mine the  state  estimates  and  associated  error  covaxiajice. 

Previous  Air  Force  Institute  of  Technology  (AFIT)  thesis  projects 
by  Robert  Mitchell  (Ref  U)  and  Robert  Mann  (Hof  5)  established  a base- 
line of  information  for  this,  problem.  In  particular,  the  thesis  by 
Hunn  extended  the  work  of  Mitchell  and  utilized  an  extended  Kalman 
filter  to  perform  state  estimation.  Mann  conducted  a covariance 
analysis  for  two  reduced  order  Kalman  filters.  As  indicated  by  the 
Mann  covariance  analysis,  the  performance  of  the  two  filters  was  not 
satisfactory.  However,  a covariance  auialysis  is  viewed  as  only  a 
first  step  in  the  evaluation  of  the  proposed  reduced  order  filters. 

The  covariance  equations  provide  filter  error  statistics  but  the 
covaxi2Lnce  analj'sis  does  not  represent  a complete  system  simulation. 

A direct  statistical  simulation  (Monte  Carlo  analysis)  is  required  to 
continue  analysis  of  the  filter  performance.  In  older  to  limit  the 
scope  of  this  study,  a Monte  Carlo  analysis  of  the  Mann  six  state 
filter  will  be  performed  since  it  is  more  likely  to  be  Implemented  if 
it  meets  performance  objectives. 

Assumptions  and  Limitations 

Since  the  system  dynamics  aare  non-linear  for  the  aircraft  to 
satellite  tracking  problem,  the  basic  Kalman  filter  cannot  be  applied 
in  this  study.  Several  non-linear  estimation  methods  axe  available 
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for  handling  problems  of  this  type.  However,  as  long  as  performance 
is  satisfactory  the  extended  Kalman  filter  will  be  used  because  the 
equations  are  simple  and  easy  to  calculate. 
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The  six  state  filter  will  be  analyzed  with  the  tracking  profile 
used  in  the  Majrin  thesis.  With  this  profile,  the  tracker  (aircraft), 
at  the  beginning  of  the  simulation,  lies  in  the  orbit  plane  of  the 
satellite.  As  the  simulation  progresses,  the  aircraft  moves  orthogonal 
to  the  orbit  plane.  The  assumption  is  made  that  this  profile  repre- 
sents a worst  case  condition.  This  is  a reasonable  assimption  because 
the  tracking  geometry  restricts  the  flow  of  information  to  some  of 
the  states  in  the  extended  Kalman  filter  and  creates  observability 
problems.  The  six  state  reduced  order  filter  will  be  tuned  using  this 
tracking  profile. 

The  remaining  assumptions  parallel  those  of  Haunt 

1.  Essentially  perfect  measurements  of  the  tracker  acceleration 

j 

with  respect  to  inertial  space  are  available  as  the  derived 
output  from  three  accelerometers, 

2.  The  tracking  system  will  provide  noise-corrupted  meaLSurements 
of  the  inertial  angular  rates  of  the  tracker,  range  and  small 
angle  deviations  between  the  boresight  and  line-of -sight 

frames.  j 

> 

3.  The  tracking  system  can  be  instantaneoxisly  corrected  by  the  j 

1 

extended  Kalman  filter  state  error  estimates.  j 


4. 

5. 


Essentially  perfect  measurements  of  the  tracker  elevation  and 


azimuth  angles  are  available  from  resolvers. 

The  tracker  y axis  will  be  inertlall5’-  stabilized  so  that  it 
always  lies  parallel  to  the  geocentric  x-y  plane. 


Research  Stratep:.v 

The  research  for  this  project  is  divided  into  the  following  five 
major  categories j 

1.  To  investigate  the  background  of  the  problem. 

2.  To  determine  the  computer  simulation  method  for  the  general 
Monte  Carlo  performance  analysis. 

3.  To  develop  a general  llonte  Carlo  Analysis  Program  (l^CAp) 
to  perform  the  system  simulation. 

4.  To  perform  the  llonte  Carlo  analysis  of  the  six  state  reduced 
order  Kalman  filter  using  HGAP. 

5.  To  analyze  the  results  of  the  computer  simulation. 

Summary 

This  study  is  concerned  with  the  development  of  a general  Monte 
Carlo  analysis  program  and  its  application  to  a problem  previously 
investigated  in  a limited  manner  throxigh  covariance  analysis.  Chapter 
II  presents  the  Kalman  filter  and  extended  Kalman  filter  equations 
used  in  this  study.  Chapter  III  briefly  describes  the  Monte  Carlo 
analysis  simulation  method  and  how  the  error  statistics  are  computed. 

A more  detalxed  description  of  HCAP  is  contained  in  the  users  guide 
(Appendix  A).  The  fourth  chapter  presents  a detailed  description  of 
the  satellite  and  aircraft  tracker  system  state  equations  for  the 
problem  investigated  with  the  general  Monte  Carlo  analysis  program. 
Chapter  V presents  the  development  of  a reduced  order  extended 
Kalman  filter  model  used  in  the  analysis.  The  results  of  the  analysis 
are  presented  and  discussed  in  Chapter  VI . The  final  chapter  presents 
the  conclusions  and  recommendations  of  this  project. 
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II,  Extended  Kalman  Filter  Equations 

Introduction 

This  chapter  presents  the  propagation  and  update  equations  for 
both  linear  and  extended  Kalman  filters,  A Kalman  filter  is,  concisely 
stated,  an  optimal  recursive  data  processing  algorithm  for  the  deter- 
mination of  the  states  or  parameters  of  a system  using  noise  corrupted 
measurements  (Ref  6:l6).  If  a physical  system  of  interest  can  be 
modeled  by  a set  of  ordinary  linear  differential  equations  and  linear 
measurements  with  system  and  measurement  noises  which  are  white  and 
Gaussian,  then  the  Kalman  filter  will  provide  the  best  estimate  of 
the  system  states.  However,  in  many  cases  of  practical  interest, 
physical  systems  must  be  represented  by  a nonlinear  set  of  differential 
equations  and/or  nonlinear  measurement  equations.  For  such  problems, 
it  is  often  convenient  to  linearize  the  system  equations  about  some 
assumed  set  of  nominal  conditions  and  use  developed  algorithms  (such 
as  the  extended  Kalman  filter  which  uses  reevaluation  of  the  nominal 
at  each  measurement  time)  for  estimation  about  these  nominal  conditions 
(Ref  7«57). 

Considering  the  approximations  necessary  aiid  the  fact  that  there 
is  no  "best"  suboptimal  filter,  the  extended  Kalman  filter  gain  and 
covariance  propagation  equations  have  the  same  form  as  the  Kalman 
filter  equations,  but  axe  linearized  about  the  current  state  estimates. 
The  linearization  is  a first  term  approximation  to  a Taylor  series 
expansion  to  the  state  estimate.  Higher  order  and  more  exact  approxi- 
mations can  be  achieved  by  using  more  terms  of  the  Taylor  series 
expansion  for  the  nonlinearities  and  by  deriving  approximate  recursive 
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relations  for  the  higher  moments  of  the  state  vector  x (Ref  8il84), 
yielding  higher  order  nonlinear  filters.  As  might  be  expected,  if 
the  system  nonlinearities  are  significant,  then  neglecting  the  higher 

i 

order  terms  will  result  in  biased  estimates.  However,  when  compared 
to  the  extended  Kalman  filter,  the  higher  order  filters  are  both  more 
complex  and  more  costly  in  texins  of  computer  implementation.  For  this 
reason,  the  extended  Kalman  filter  is  often  considered  first  in  non- 
linear estimation  problems. 

Notation 

This  study  has  adopted  the  notation  presented  by  Wrigley,  Hollister, 
and  Denhard  (Ref  9*20-23).  A vector  (represented  by  a letter  with  an 
underbar,  x)  is  considered  to  be  a geometric  entity  in  real,  three 
dimensional  space.  The  vector  represents  some  physical  quantity  which 

i 

has  both  magnitude  and  direction.  When  the  physical  quantity  is  « 
measured  with  respect  to  a coordinate  frame,  the  vector  is  said  to  be 
coordinatized  in  that  frame.  The  three  numbers  associated  with  the 
mathematical  vector  cire  the  components  of  the  physical  vector  relative 
to  the  specified  coordinate  frame.  As  an  example,  if  x is  coordinatized 
in  the  "i"  ilame,  the  vector  would  be  denoted  by  x^,  a tlaee- tuple  of 
numbers.  Another  vector  of  interest  is  the  physical  angular-velocity 
vector,  generally  denoted  by  a w.  The  angulair  velocity  vector  will 
have  two  subscripts,  as  an  example  The  subscripts  indicate  that 

the  angular  velocity,  w,  is  of  frame  b with  respect  to  frame  n, 
coordinatized  In  the  1 frame. 

1 

j 
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A 3x3  direction  cosine  matrix  is  used  in  this  study  to  transform 
the  components  of  a vector  in  one  frame  to  those  in  another.  The 
direction  cosine  matrix  is  represented  by  a capital  C and  an  undertar. 
Associated  with  the  letter  C will  he  a subscript  for  the  frame  from 
which  the  transformation  is  made,  and  a superscript  for  the  new  coordi- 
nate frame.  As  an  example,  would  transform  the  vector  x from  the 
b frame  to  the  n frame,  as  x*^  = x^, 

V/here  it  is  necessary  to  address  individual  components  of  a vector 
coordinatized  in  a specific  frame,  the  vector  will  be  specified  and 
subscripts  used  to  indicate  individual  components.  For  example: 


w = 


w. 


w; 


i i iT 

Wx  wy  w^ 


where  T denotes  the  transpose  operator. 


(1) 


Definitions 

Listed  below  are  some  of  the  definitions  used  in  this  chapter: 

x(t^)  = system  state  at  time  tj^  (n-vector) 

x(t^)  = filter  estimate  prior  to  incorporating  a measurement  at 
time  tj^  (n-vector) 

x('t^)  = filter  estimate  after  incorporating  a measurement  at  time 
tj^  (n-vector) 

^(tj^)  = m vector  of  measurements  at  time 

w(-t^)  = system  dynamics  white  Gaussian  noise  s-vector,  independent 
of  x(tQ),  where  Cw('t)]  “ 0»  and  iw(t)w(T)'^  “ ^(t)d(t-T), 
[[w(t)  is  assumed  to  be  zero  mean,  Gaussian,  and  white 
(uncorrelated  in  time)^  with  ^(t)  an  sxs  positive  semi- 
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definite  symmetric  matrix  that  is,  in  general,  piecewise 
continuous  in  t. 

v(tj^)  = zero  mean,  white  Gaussian,  measurement  noise  sequence 
independent  of  w(t)  and  xCtg)  for  all  time  (m-vector). 

The  statistics  of  v(tj^)  are  C v(tj^)]  = 0,  and 
iv(tj^)v(tp]  = R(t^)  = tj 

0 otherwise 

= state  transition  matrix  from  time  to  time  tj^ 

(nxn  matrix) 

P(t^)  = filter  computer  covariance  matrix  of  state  x(tj^),  also  of 
the  error  in  the  estimate  of  prior  to  incorporating 

a measurement  at  time  (nxn  symmetric  matrix) 

P(t^)  = filter  computer  covariance  matrix  of  state  x(tj^)».also  of 
the  error  in  the  estimate  of  x('ti)»  after  incorporating  a 
measurement  at  time  tj^  (nxn  symmetric  matrix) 

F(t)  = system  dynamics  matrix  (nxn)  or  the  matrix  of  partial 

derivatives  of  iCx('t)»tl  with  respect  to  x for  the  extended 
Kalmcin  filter 

G(t)  = system  input  matrix  (nxs) 

H(tj^)  = system  observation  matrix  at  time  tj^  (mxn)  or  the  matrix 
of  partial  derivatives  of  l£x(t),t!l  with  respect  to  x 
for  the  extended  Kalman  filter 

^(■tj^)  “ positive  definite  measurement  noise  covariance  matrix  (mxm) 

K(tj^)  = Kalman  gain  matrix  (nxm)  defined  at  time  t^^ 
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Linear  Kalman  Filter  Formulation  \ 

The  linear  Kalraan  filter  formulation  presented  in  this  section  is 
for  a continuous  time  system  model  with  discrete  time  measurement 
updates.  Assume  that  the  system  modeling  has  been  completed  and  that 
the  state  vector  x(t)  satisfies  the  vector  stochastic  differential 
equation! 

iCt)  = F(t)  x(t)  + G(t)  w(t)  (2) 

The  state  equation  is  propagated  forward  in  time  from  the  Initial 
condition  x('to)*  Since  the  exact  Initial  condition  may  not  be  known, 
it  is  modeled  as  being  a Gaussian  random  variable  with  mean  and 
covariance  P : 


0^ 

1 1 

n 

(3) 

(^) 

The  initial  covariance  matrix  may  be  positive  semidefinite, 
admitting  exact  knowledge  of  the  initial  conditions  of  some  of  the 
states  or  linear  combinations  thereof.  It  can  be  shown  (Ref  10il57- 
163)  that,  under  the  assujnption  that  x(to)  is  either  deterministic 
or  a Gaussian  random  variable,  the  solution  x('t)  linear  stochastic 
differential  equations  such  as  Equation  (2)  is  a Gauss-Harkov  process, 
l.e,  the  conditional  density  of  x at  time  tp  based  upon  all  realiza- 
tions of  X through  time  tj,  is  both  Gaussian  and  completely  determined 
by  the  process  value  at  ti-i.  Because  the  conditional  density  is 
Gaussian,  it  is  completely  specified  by  its  mean  and  covariance 
(Ref  7»92). 
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Measurements  are  available  at  discrete  time  points  and  eire 
assumed  to  be  of  the  form  of  a linear  combination  of  the  states  and 
corrupted  by  a white  Gaussian  sequence  (Jief  }s2)t 

' H(tj^)  x(ti)  + v(t^)  (5) 

The  state  estimate  propagates  between  measurements  (from  time 
ti_i  to  time  tj^)  according  to; 

i(\)  i(’tt-i)  (6) 

and  the  covariance  propagates  according  to; 

P(ti)  P(tt-i) 

+ S ^ t(ti,T)  g(t)  'i'T  (7) 

^i-l 

At  neasiurement  time  t^^,  the  estimate  is  updated  according  to 
(Ref  10; 233); 

i(tt)  = x(ti)  + K(t^)  [1^  - H(t^)  |(tp]  (8) 

P(tt)  = P(ti)  - K(t^)  H(ti)  P(tp  (9) 

where 

U\)  - P(t:)  H(tj^)’’  [H(t^)  P(tp  H(ti)^  + R(t^)r^  (10) 

where  [ indicates  the  inverse  of  the  bracketted  matrix  and  li  is 
the  realised  .'aluc  of  the  measurement  z(tj^)  at  time  tj^. 

Under  the  assumption  that  the  adequate  system  model  is  linear, 
and  that  the  dynamic  driving  and  measurement  noises  are  Gaussian  and 
white,  the  Kalman  filter  provides  the  optimal  estimate  x(ti)  of  the 
state  of  the  system  (Ref  10;66,214),  relative  to  many  optimality 
criteria  with  these  assumptions  x(t^)  is  the  mean,  mode,  and  median 
of  the  conditional  density  of  x(ii)»  conditioned  on  the  entire  measure- 
ment history  through  time  tj^.  The  covariance  of  the  error  committed 

13 


by  usin^j  x(t^)  as  the  estimate  of  the  state  at  time  tj^  is  denoted  by 
P(t^),  It  should  be  noted  that  for  a linear  estimation  problem,  the 
covariance  propagation  C^iuations  (?,9)1»  while  depending  on  the 
sequence  of  H(tj^)  and  R(tj^),  is  independent  of  the  time  history  of 
measurements  X2*  •••)•  This  will  no  longer  be  the  case  in  the 

extended  Kalman  filter  formulation. 

The  assumption  that  the  system  can  be  modeled  as  being  driven  by 
white  Gaussian  noise  is  often  well  founded  on  two  accounts.  First, 
it  has  been  found  from  practical  experience  that  the  Gaussian  distri- 
bution provides  a reasonable  approximation  to  observed  random  behavior 
in  certain  physical  systems  (Ref  7»92),  Secondly,  the  central  limit 
theorem  (Ref  7*96)  states  that  if  the  random  phenomenon  that  we  observe 
at  the  macroscopic  level,  is  due  to  the  superposition  of  a large  number 
of  independent  random  variables  on  the  microscopic  level,  then  the 
macroscopic  phenomenon  can  be  adequately  modeled  as  a Gaussian  random 
vaariable  (Ref  10*40), 

Extended  Kalman  Filter  Formulation  (Ref  3*179-189) 

The  extended  Kalman  filter  formulation  is  commonly  used  in  esti-  || 

mation  problems  in  which  the  adequate  stato  and/or  measurement  equations 
axe  nonlinecor  rather  than  linear.  Consider,  as  before,  a system  that 
is  continuous  in  time  with  measurements  at  discrete  sampling  times. 

Assume  that  the  system  state  satisfies  the  following  nonlinear  vector 
stochastic  differential  equation* 

i(t)  - fCx(t),  t]  + G(t)  w(t)  (11) 
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where  is  a nonlinear  fxmction  of  the  states  and  time  [^in 

general  f could  also  be  a function  of  deterministic  control  inputs 
ii(t)],  and  where  the  vector  w(t)  of  zero  mean  white  Gaussian  driving 
noises  and  covariance  ^(t)  enters  in  a linear  additive  manner.  The 
initial  condition  of  x(tQ)  is  modeled  as  a Gaussian  random  variable 
with  mean  and  covariance  P^.  Noise  corrupted  vector  measurements 
of  a nonlinear  function  of  the  states  and  time  are  available  at 
discrete  times  t^  asi 

z(ti)  = + v(tj^)  (12) 

where  v(tj^)  is  a zero  mean  vrhite  Gaussian  sequence  with  covariance 
kernel  (Ref  3*120) • 


Cx('ti) 


(13) 


R(tj^)  for  i = j 

0 otherwise 

It  is  assumed  that  the  processes  w and  v are  independent  of  each  other. 
The  filter  propaigation  eqiaations  are: 

^ i(Vti)  =£i(  ■!/%),  t]  (14) 

i(ti/tj^)  = x(t^)  (15) 

where  the  notation  x(t/tj^)  means  the  optimal  estimate  of  the  state, 

X,  at  time,  t,  given  the  updated  estimates  up  to  and  including  time 
tj^.  In  addition,  (the  covariance  is  propagated  approximately  by): 
P(t/tj^)  = F[t;x(t/tj^)]  P(t/tj^)  + P(t/tj^)  |[t;x(t/ti)]^ 

+ G(t)  ^(t)  G(t)'^  (16) 

P(ti/ti)  - P(t^)  (17) 

where  F is  the  matrix  of  partial  derivatives  of  f with  respect  to  x, 
evaluated  along  the  trajectory  [ which  is  propagated  by  means  of  Equation 
(14)]: 
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The  measurenent  update  is  given  byj 


li('ti)  “ P(ti)  P(t7)  i£tj^;x(tp] 


+ H(ti)r^  (19) 

- i£x(ti),tj^]} 

(20) 

£(ti)  = P(%)  - K(t^)  KLt;x(tp]?(t7) 

(21) 

where 

(22) 

P(ti)  = P(ti/t^.i) 

(23) 

r T 

>£t,x(q)]  - ^ 

■ (24) 

The  notation  t^  implies  the  value  of  the  quantity  at  the  instant  prior 
to  the  update  at  time  tj^,  and  t^  is  the  value  of  the  quantity  just  after 
the  update.  The  notation  used  thus  fair,  except  as  noted  for  the  F 
matrix,  is  developed  and  used  in  (Ref  3)« 

Unlike  the  conventional  Kalman  filter,  the  extended  Kalman  filter 
||  gain  and  estimation  error  covariance  matrices  depend  on  the  time  history 

of  x(t/tj^).  Since  a covariance  analysis  could  only  use  an  ^nominal^^^ 
as  an  approximation  to  x(t/tj^)  as  actmlly  used  by  the  filter  in  its 
online  operation  ( the  linearized  Kalman  filter) , the  extended  Kalman 
filter  performance  should  be  verified  by  a Ilonte  Carlo  simulation. 

In  general,  this  form  of  analysis  is  both  more  time  consuming  and  more 
costly  in  terms  of  computer  usage  than  is  a covariance  analysis. 
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Ill,  Monte  Carlo  Simulation  and  Analysis  of  a Sub-Optimal 

Kalman  Filter 


Introduction 

As  mentioned  in  Chapter  I,  a Monte  Carlo  analysis  is  a computer 
simulation  technique  that  c' •<  be  used  to  develop  the  required 
(extended)  Kalman  filter  w .ing,  portray  filter  performance  capabilities, 
and  study  the  effects  of  parameter  variations  in  the  true  system  model. 
The  ilonte  Carlo  simulation  uses  random  number  generators  and  shaping 
filters  which  generate  random  errors  do  noise  corrupt  both  system 
state  dynamics  and  measurements.  The  noisy  system  measurements  are 
then  processed  by  the  suboptimal  (extended)  Kalman  filter  to  generate 
the  filter  estimates.  If  many  of  these  simulations  are  conducted, 
then  one  can  compute  the  statistics  of  the  error  between  the  filter 
estimate  and  the  truth  model  for  quantities  of  interest. 

Performance  Analysis 

Figure  1 depicts  schematically  a means  of  conducting  the  perform- 
ance analysis  of  a given  (extended)  Kalman  filter  design.  The  truth 
model  by  definition  is  the  best,  most  complete  mathematical  model  that 
can  be  developed  to  describe  the  system  under  study.  Such  a truth  model 
is  the  product  of  extensive  study  and  data  analysis  of  the  system. 

As  noted  in  the  figure,  the  truth  model  is  an  n-fState  model,  linear 
or  nonlinear,  driven  by  noise  w-t(t)  (assumed  white  and  Gaussian),  that 
generates  the  true  state  values  x^(t).  It  is  assumed  that  the  true 
values  for  the  critical  quantities  of  interest  eore  related  to  the  true 
state  values  by  the  vector  function  (a  p vector  - generally  less 
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than  n^)  i 

(25) 

Often,  c^Cx^(t),t]  is  a linear  function  of  x^,  as 

1^(^)  “ £(t)x(t)  (26) 

In  fact,  in  many  instances,  the  quantities  of  interest  are  simply  a 
subset  of  the  truth  model  states;  letting  the  components  of  x(t)  be 
ordered  so  that  the  first  p states  are  the  quantities  of  interest, 
this  yields  the  particular  formi 

= Ci|^  st^t)  (27) 

Thus  we  have  access  to  the  "real  world"  values  for  comparison  with  the 
filter  estimates. 

At  discrete  times  {t^)  the  m-dimensional  measurements* 

are  presented  to  the  Kalman  filter.  The  Kalman  filter  algorithm 
processes  the  measurements  and  produces  the  state  estimates  x(t). 

These  state  estimates  are  related  to  the  estimates  of  critical  quanti- 
ties of  Interest  by  the  vector  function  c (a  p vector): 

i(t)  =c[5(t),t]  (29) 

Thus  the  true  error  committed  by  the  Kalman  filter  (£^)  at  time  t^, 
before  and  after  measurement  incorporation,  is* 

(30) 

£t(tt)  - i(tt)  - (31) 

If  a feedback  control  system  with  impulsive  correction  is  assumed, 
then  a third  error  is  of  interest  as  well* 

- Htf)  - it(ti°)  (32) 

where  the  superscript  c denotes  after  the  control  is  applied. 


'i 

i 


The  objective  of  the  performance  analysis  is  to  characterize  the 
error  process  statistically.  This  is  accomplished  in  the  ilonte  Carlo 
simulation  by  generating  many  samples  of  the  error  process  and  then 
computing  the  sample  statistics  directly.  If  enough  samples  are  gene- 
rated, then  the  sample  statistics  should  approximate  the  process 
statistics  very  well  (Ref  The  sample  statistics  con.puted 

at  each  point  in  time  are  the  calculated  mean  error  (je),  calculated 
standard  deviation  of  the  error  (s^),  ?nd  the  ensemble  average  of  the 
filter  co\^riance  diagonal  terms  . The  calculations  are  made 

over  the  ensemble  of  runs  (n)  for  each  time  point  (tj)» 


i(tj)  = N ^2^[xj^(t^)  - iti^^j^^ 

Sg(t)  = [c^Z^Cx^Ctp  - Xti(^j)]^  ■ N (3^) 

^(tj) -i  (35) 

These  equations  can  be  found  in  any  good  statistics  text  (see  for 
example  Ref  lli26,40).  It  is  important  to  note  the  use  of  l/(K-l) 
in  the  expression  for  the  standard  deviation,  Sg,  instead  of  l/ll; 
this  results  in  an  unbiased  estimator  of  variance  before  the  [ ^T/2 
is  taken  and  is  customarily  used  for  sampled  data  (for  small  or 
moderate  N), 

One  useful  output  of  the  Monte  Carlo  analysis  is  a piot  of  the 

filter  estimate  of  the  standard  deviation,  t . ) . along  with  the 

«3 

corresponding  computed  standard  deviation,  £e(t.)  for  all  t^  or 
Interest.  The  Kalman  filter  is  timed  by  minimizing  the  computed  error 
standard  deviations.  As  a rule  of  thumb,  good  tuning  is  achieved 


I| 


:| 

I 


20 


i 
1 

when  the  filter  estimated  standard  deviations  are  approximately  equal 

to  the  computed  standard  deviations.  If  the  problem  under  investlga-  j 

tion  is  sensitive  to  small  changes  in  the  filter  tuning  parameters 
£p(tp),  and  time  iiistories  of  i^(t)  and  H(t)]  then  the  investigator 

must)  ! 

1,  Conduct  the  Monte  Carlo  simulation 

2,  Examine  the  standard  deviation  plots  i 

3,  Determine  tuning  parameter(s)  to  be  changed  I 

4,  Repeat  1 through  3 until  adequate  tuning  is  achieved 
Thus  it  can  be  readily  visualized  that  precise  timing  of  the  Kalman 
filter  by  Monte  Carlo  techniques  is  very  costly  in  terms  of  effort  and 
time. 

A second  useful  output  of  the  Monte  Carlo  analysis  is  a plot  of 
||  the  mean  en'or  versus  time.  This  plot  is  usually  obtained  with  one 

L standard  deviation  bounds  plotted  above  and  below  the  mean  error, 

I This  plot  is  useful  to  determine  if  the  suboptimal  Kalman  filter 

It  provides  the  desired  accuracy.  If  the  filter  starts  from  initial 

conditions  other  than  the  truth  model  initial  conditions,  then  the  ' 

filter  performance  to  non-zero  initial  error  can  be  observed.  Bias  I 

errors,  a problem  of  significance  to  many  extended  Kalman  filter  j 

designs,  also  become  readily  apparent  from  the  plots  of  this  form. 

Digital  Simulation  of  the  Truth  Model  and  Filter  Model 

In  order  to  simulate  the  system  on  a digital  computer,  the  problem 
must  be  divided  into  segments  which,  when  processed  sequentially  over 

a given  time  period,  will  provide  a realistic  portrayal  of  filter  use  ^ 
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in  actual  on-line  implementation.  A convenient  time  period  5s  the  i 

sample  period  between  times  at  which  measurements  become  available.  | 

Then  the  entire  procedure  is  iterated  over  some  desired  time  interval  j 

of  interest.  The  segments  arej 

1.  Time  propagation  of  system  state  eq,uations 

2.  Time  propagation  of  filter  state  equations 

3.  Tine  propagation  of  filter  covariance  equations 

4.  Update  of  filter  covariance  matrix 

5.  Generation  of  measurements  for  the  filter 

6.  Update  of  filter  states 

7.  Application  of  impulsive  feedback  to  "update"  system  states 

8.  Reset  the  filter  states  after  the  feedback 

The  first  segment  is  simulated  by  integrating  the  n^-dimensional 
system  (truth  model)  stochastic  differential  cqiiationj 

= f^[x^(t),t]  + G^(t)w^(t)  (36) 

from  time  tj^  to  time  of  the  next  measurement  time  tj^+^.  The  integra- 
tion is  typically  performed  by.  either  auler  or  Runga-Kutta  numerical 
integration  methods  with  a specified  integration  step  size.  Integrating 
the  deterministic  part  of  the  differential  equation  is  readily  apparent,  • 

however,  the  additive  white  Gaussian  noise  w.j^  is  not.  The  following 
derivation  shows  how  the  solution  to  the  stochastic  differential 
equation  can  be  approximated. 

To  provide  insight  into  the  approximate  solution  form,  consider 
firsw  a linear  stochastic  differential  equation  of  the  formi 

i(t)  = F(t)  x(t)  + G(t)  w(t)  (37) 
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where  w(t)  is  a zero  mean  white  Gaussian  noise  with  i:(^w(t)w(t+t)'^]  = 

^(t)6(t).  An  equivalent  discrete  time  system  equation  would  bej 

^][  + l 

x(  ti+i)  = i(  ) x(  t t'i+ifti)  g(t)w(t)  dr 

h (38) 

= 4(Vi>^l^  * ^(^i)  (39) 

where  0(*,*)  satisfies  both  of  the  following  relationships* 

^ ^("titj^)  = E.(t)  2(t,tj^)  for  t in  the  interval 

(ti.Vi)  (^0) 

$(ti.t^)=i  {ki) 

For  this  discrete  time  system,  i!:jj^(‘tj^)  is  a zero  mean  white  Gaussian 
discrete-time  noise  with  strength* 

Cwd(\)iid(^i)^l  = / ^ ^ ^(ti+i,T)  G(t)  RCt)  G(t)’' 

ti  (42) 

^(tj^+^.r)  dT 

= 2d(ti)  (43) 

2[wd('’^i)iJd(^j^^^  " ^ ^i  (^) 

How  if  the  integration  step  size  is  small  relative  to  system  character- 
istic response  tine  and  constant  (At  sec  in  duration)  then  the  following 
approximations  can  be  made: 

i(ti  i.t^)  = I + F(t^)  At  (45) 

i G a dT  = [G(t^)  a(ti)  G(t^)‘^]  At  (46) 

^1 

2 

where  all  terms  of  order  At  or  higher  have  been  neglected.  Thus  one 
approximate  simulation  of  the  continuous  tine  system  would  be* 

i(ti+i)  “ [l  + F(tj^)  4t]  xit^)  + ^i\)  (47) 

= x('ti)  + ll(ii)  ^i\)  + Wd(\)  (^) 


Khere  is  described  by: 

a(\)  G(t^)'^l  At  (49) 

Thus  the  differential  equation  is  integrated  usin,:5  Equation  (46) 
and  noise  generators  utilized  for  the  w^(tj^)  terra.  Higher  order  inte- 
gration techniques  can  be  applied  to  the  nonhomogeneous  portion  of  the 
differential  equation,  with  the  same  w^(t^)  as  above  used  as  the  forcing 
terms.  It  should  be  noted  that  this  technique  is  directly  extendable 
to  nonlinear  equations  of  the  form: 

iCb)  = l[x(t),t]  + G(t)  w(t)  (50) 

Segment  two  is  performed  by  calculating  the  filter  equation  (6). 

The  filter  covariance  propagation  Equation  (?)  is  calculated  to 
accomplish  segment  three.  Segment  foiur  is  perfonaed  by  solving 
Equations  (lO  and  9).  Segment  five  uses  Equation  (5)  to  generate  the 
m-dimerisional  measurements.  Noise  generators  are  utilized  to  provide 
the  v^(t^)  term.  The  filter  update  Equation  (8)  is  calculated  to 
perform  segment  six. 

A Kalman  filter  will  often  be  implemented  in  an  indirect  feedback 
configxuration,  in  which  filter  error  state  estimates  are  fed  back  to 
the  actual  system  to  correct  the  actual  error.  If  the  system  correc- 
tion can  be  performed  rapidly  enough  compared  to  the  filter  update 
period,  then  the  feedback  can  be  simulated  as  an  impulsive  change. 

After  application  of  the  correction,  assaming  linear  feedback,  the 
truth  model  state  process  becomes  (Ref  3*6-76): 
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The  filter  should  be  "told"  that  this  feedback  to  the  system  has 
occurred,  so  its  state  estimate  is  modified  asi 

x(ti®)  = «(t^)  - D(t^)  x(t^)  (52) 

= [i  - D(tj^)]  x(ti)  (53) 

Note  that  if  the  filter  is  estimating  error  states  then  the  U(tj^) 
matrix  will  typically  be  an  identity  matrix.  Segments  seven  and 
eight  are  performed  by  calculating  the  results  of  Equations  (51) 
and  (52). 

Each  of  these  segments  are  performed  sequentially  at  the  measure- 
ment update  time  and  the  simulation  time,  system  state,  filter  state, 
and  filter  covariance  values  saved  on  a data  tape.  This  process  con- 
tinues until  the  simulation  final  time  is  reached,  thus  completing 
one  run  of  a Monte  Carlo  simulation.  After  many  runs  of  the  iionte 
Ceirlo  simulation  have  been  performed,  then  the  data  tape  can  be  read, 
the  sample  statistics  calculated  using  Equations  (33”35)  and  computer 
plots  created. 
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IV.  Satellite  and  A^ircrr^ft  Tracker  System  State  Equations  i 
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I 

1 

Introduction  j 

In  this  chapter,  the  satellite  (target)  and  aircraft  tracker  sys- 
tem state  and  measurement  equations  are  developed  for  the  system 
(truth)  model.  Since  this  study  is  a follow-on  to  the  work  of  Mitchell 
(Ref  4)  and  Mann  (Ref  5)#  the  presentation  in  this  chapter  will  paral- 
lel their  development,  i 

System  Description 

Physically,  the  problem  consists  of  an  orbiting  satellite  (target) 
and  a moving  aircraft  which  is  equipped  with  some  type  of  tracking 
device.  The  tracking  device  is  typically  a radar  or  laser.  Because 

the  aircraft  ajid  satellite  are  both  moving  in  this  problem,  the  i 

tracking  system  must  be  capable  of  following  the  satellite  motion.  | 

The  tracking  system  is  equipped  with  three  rate  gyros  which  provide  j 

I 

measurements  of  the  three  components  of  the  tracker  Inertial  angular  | 

velocity.  Because  of  tracking  system  dynamics  and  errors  in  the  j* 

I 

tracking  system,  the  satellite  will  not  lie  exactly  along  the  bore-  ' 

sight  of  the  tracking  device.  In  this  problem,  a tracker  control  sys-  i 

tem  is  available  to  correct  the  estimated  tracker  error  angles.  These  j 

estimates  of  the  error  angles  will  be  provided  hy  an  extended  Kalman 
filter.  Three  accelerometers,  mounted  on  the  aircraft,  provide 

measurements  of  the  specific  force  of  the  tracker.  These  specific  ; 

i 

force  measurements  are  used  hy  the  aircraft  Inertial  Navigation  Sys- 
tem (ins)  to  determine  the  tracker  Inertial  position  amd  velocity.  ; | 


r 


In  this  study,  It  Is  assuned  that  uncorrupted  (treated  as  perfect) 
tracker  Inertial  position  and  velocity  information  is  available  from 
the  aircraft  INS.  In  addition  to  the  tracker  angular  velocity,  the 
tracking  system  provides  measurements  of  range  to  the  target  and  the 
true  line  of  sight.  The  measurements  are  all  assumed  to  be  imperfect. 
The  nature  of  the  noise  corrupting  the  measurements  will  be  discussed 
later. 

Satellite  State  So nations 

The  system  model  used  for  the  satellite  dynamics  will  be  presented 
in  this  section.  The  satellite  state  equations  will  not  be  formally 
derived;  the  approach  taken  is  straightforward  and  the  derivations  of 
the  models  may  be  found  by  the  interested  reader  in  any  good  astro- 
dynamics  text  (see  for  example  Ref  12),  As  noted  in  Mann  (Ref  5 ‘8), 
due  to  limitations  in  the  computer  simulation,  the  solar  pressure 
perturbative  acceleration  is  deleted  from  Mitchell's  system  model.  To 
compensate  for  this  additional  unmodeled  effect,  the  strengths  of  the 
driving  noise  sources  on  the  satellite  equations  are  Increased  by  an 
appropriate  amount. 

The  target  state  equations  are  expressed  in  the  geocentric 
equatorial  nonrotating  coordinate  system  with  the  x-axis  lying  along 
the  line  of  the  mean  vernal  equinox  (Figure  2),  This  coordinate  sys- 
tem will  be  considered  to  be  aji  Inertial  ("i”)  frame  for  this  appli- 
cation. The  state  equations  describing  the  satellite's  motion  are; 


where 
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(56) 

(57) 

(58) 

(59) 

Xg,  represent  the  target  inertial  position  components 
along  the  x,  y,  z axes  respectively 
x^,  xg  represent  the  target  Inertial  velocity 
a is  the  earth’s  gravitational  acceleration  vector 
a^ij  is  the  Ixmar  gravitational  perturbation  vector 
a^  is  the  solar  gravitational  perturbation  vector 
^ is  the  atmospheric  drag  acceleration  vector 
Wi,  Wg,  w^  are  zero-mean  independent  white,  Gaussian  noises 
included  to  account  for  unmodelcd  effects.  These 
effects  include:  solar  pressure  perturbations, 
higher  order  gravitational  terns,  Jind  uncertain- 
ties in  the  models  used  in  this  study,  such  as 
deviations  in  the  atmospheric  density. 

In  order  to  determine  the  strengths  of  the  white  noises,  consider 
the  following  calculations.  For  a relatively  small  satellite  - in  a 
200  Km  circular  near  polar  orbit  - with  a solar  pressure  coefficient 
equal  to  that  of  a vehicle  with  a projected  surface  towards  the  sun  of 
10  m^  and  a ballistic  coefficient  of  0.15»  the  terms  of  Xji^  have  deter- 
ministic values  of: 

& “ -7.55  m/sec  - acceleration  due  to  full  gravity 
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ajjjjj  “ +5.0x10  m/sec  - lunar  perturbative  acceleration 
~9  2 

agjj  =■  +2,0x10  ra/sec  - solar  gravitational  perturbative 

acceleration 

—9  2 

^■dx  * "9.0x10  m/sec  - drag  acceleration 

These  values  have  been  determined  using  the  models  proposed  later  in 

this  section  with  the  satellite  at  north  latitude  and  the  sun  and 

moon  positioned  for  worst  case  effects.  The  unmodoled  solar  pressure 

perturbation  on  the  satellite  imler  these  assumptions  vrould  be 
-9/2 

-2.0x10  m/sec  . Because  of  the  aforementioned  criteria  of  modeling 

all  perturbative  accelerations  of  magnitude  greater  than  10  m/sec  , 

a reasonable  value  for  the  distribution  standard  deviation  of  w^  due 

-9  / 2 

to  modeling  uncertainty  and  higher  order  effects  is  1x10  m/sec  . 

Taking  into  account  the  unmodeled  effects  due  to  the  solar  pressure 

pertxirbations,  w^^,  Wg,  and  w^  are  modeled  as  zero  mean  independent 

■“9 

white  Gaussian  noises  with  distribution  one  sigma  values  of  3^10  ^ m/ 
2 

sec  . 

Gravitational  Field  liodeling 

Modeling  of  the  earth's  gravitational  field  is  accomplished  in  a 
geocentric,  equatorial,  rotating  coordinate  frame.  The  relationship 
between  this  "r"  frame  (x^,  y^,  z^)'^  and  the  inertial  "i"  (x^,  y^, 
coordinate  frame  used  in  the  previous  section  is  shown  in  Figure  2, 
The  transformation  matrix  from  the  rotating  (r)  to  the  Inertial  (i) 
frame  is  defined  asi 
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sin0 
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where 


and 


9 = 9_  + w,  t 
g le 


(61) 


9 = local  sidereal  tiine  at  t=0 

g 

Wig  = earth's  angular  rotation  rate  (7.292x10  ^ rad/sec) 
t = time 

The  potential  model  that  was  chosen  includes  tesseral,  zonal, 
and  sectorial  hai-monics  up  to  and  including  (6,6)  (Ref  12il73~l80). 

The  gravitational  potential  U(x^,  y^,  z^)  in  the  "r"  frame  is  (Ref  12: 

175) « 

u . i|^  1 1(2  I’l"'  ^ cc^_,  cos(rt.s,) 


k=2  m=0  “ ^ 


where  the  terms  in  Equation  (62)  are  defined  as  follows: 
,(ni) 


(62) 


Pj^  ''(sin0)  are  Legendre  functions: 

p[~^(sin!i()  - (1  - sin^)’"^'^^  tPjcinJ^)} 

d(sin^) 

and  Pj^(sin^)  is  the  Legendre  polynomial  with  argument  sinp^. 


m 


p/i 

is  the  mass  of  the  earth  (5.983  x 10  Kg) 


2 1^ 

kg  is  the  gravitational  constant  for  the  earth  (kg  = 3.986x10 
rn^/sec^) 

r is  the  radial  distance  from  the  earth's  center  to  the  satellite 
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^ is  the  geocentric  declination  angle  of  the  satellite 
Xg  is  the  longitude  of  the  satellite  with  respect  to  the 
prine  jiieridlan. 

^k,m  ^k,ni  harmonic  coefficients  for  the  gravita- 
tional potential  such  that  and  Sj^  q = 0 and 

the  coefficients  are  the  zonal  harmonics. 

^k,n  tenaed  tesseral  harmonics  if  m k,  m 0 

and  sectorial  harmonics  if  m = k (Ref  4j115-117). 

The  components  of  the  gravitational  acceleration  vector  along 
the  X , y , z axes  - - can  no>?  he  determined  hy 
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and  the  gravitational  acceleration  vector  in  the  inertial  frame  can 

he  determined  from 
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Lunar  and  Solar  Perturbative  Accelerations 

The  perturbative  accelerations  on  the  satellite  due  to  the  lunar 
and  solar  gravitational  fields  is  presented  in  this  section.  Since 
the  time  elapsed  during  a complete  tracking  pass  is  small  vhen  compared 
to  the  inertial  dynamics  of  the  moon  and  sun,  they  are  considered 
stationary  for  the  purposes  of  this  study. 
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In  a similar  manner,  the  perturbative  acceleration  due  to  the 

sun's  gravitational  field  is 
,i 
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The  position  vector  of  the  sun  in  inertial  frame  is  defined  as 
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eind  the  perturbative  acceleration  on  the  vehicle  due  to  the  sun's 
gravitational  field  is 


(Sr)  = 


(70) 


where 


20  3 2 

= gravitational  parameter  of  the  sun  (1.32?  x 10  m /sec  ) 
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Acceleration  Due  to  Drag 


Satellite  acceleration  due  to  drag  is  modeled  as  a fvinction  of  the 
height  above  the  earth's  siarface,  the  velocity  of  the  satellite  rela- 
tive to  the  rotating  atmosphere  and  the  vehicle  ballistic  coefficient: 


where 


(a^)^  ■=•  drag  acceleration  vector  in  inertial  frame 

^5  " '^ie^l 
^6 

B = ballistic  coefficient  of  the  satellite 

—All 

k •=  atmospheric  density,  modeled  as  X.  = X^e 

= mean  sea  level  atmospheric  density  (1.376229  x 10  ^ g/m^) 

0 = altitude  atmospheric  density  decay  rate  (1.395  x 10  Vni) 

/ 2 2 2.1/2 

h = (xj^  + Xg  + x^;  - *=  height  above  mean  earth  radius 

Rq  “ mean  earth  radius  (6,37817  x 10^  m) 

Wj^g  “ angular  rotation  rate  of  the  earth 
While  the  ballistic  coefficient  of  the  vehicle  is  generally  not 
known,  it  is  known  that  for  a nonthrusting  vehicle  it  will  not  change 
significantly  during  the  time  of  a tracking  pass  (an  attitude  maneuver 
could  affect  it  by  changing  the  surface  area  along  the  velocity  vector). 
Hence,  it  is  assu.aed  in  this  study  that  the  ballistic  coefficient  can 
be  adequately  modeled  as  a random  constant  or  bias  (a  random  variable 
that  has  lOQ^o  correlation  in  time) 

B “ 0 (72) 

with  an  initial  condition  as  a Gaussian  random  variable. 

Tracker  State  Equations 

While  the  target  state  equations  are  straightforwaird  and  represent 
a commonly  used  model  for  satellite  dynamics,  the  tracker  state 
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dynamics  and  measurement  equations  are  very  dependent  upon  the  modeling 
assumptions  made  in  this  study.  Therefore,  a full  development  of  the 
tracker  dynamic  state  equations  and  then  the  tracking  system  measurement 
equations  will  be  given  in  this  section. 


The  geometry  of  the  tracker  is  shown  in  Figure  3> 


It  is  assumed  that  the  tracker  base  is  inertially  stabilized  such  ] 

1 

that  the  tracker  elevation  axis,  y , always  lies  in  the  inertial  X-Y  | 

I 

plane.  The  axis  is  along  the  boresight  of  the  tracker  and  the  z^  ; 

axis  completes  a right-hand  orthogonal  system.  Assuming  that  the  , | 

I 

tracker  base  is  inertially  stabilized  as  above,  the  two  Euler  angle  j 

rotations  needed  to  go  from  the  inertial  to  the  tracker  frame  are  | 

dependent  only  upon  the  relative  position  vector  from  the  tracker  to  I 

the  target,  expressed  in  inertial  coordinates.  The  first  Euler  angle  | 

i 

rotation  is  ty  an  angle  0 about  the  inertial  z axis  as  shown  in  Figure  4.  | 

The  superscript  "a”  indicates  an  intermediate  frame  and  the  transforma-  | 

tion  matrix  is 
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3) 


Figure  4.  First  iiuler  Angle  Rotation 
Because  of  the  constraint  on  the  axis,  the  next  Euler  angle  rota- 
tion is  by  an  angle  ^ about  the  axis  - leaving  the  y^  axis  in  the 
inertial  xy  plane,  A view  of  this  rotation  is  shown  in  Figure  5,  The 

transformation  matrix  between  the  intermediate  "a"  fraae  and  the  tracker 

"t"  frame  is 

cosp!  0 -sinJZJ 

0 10 
sinj^  0 cos^ 

Therefore,  the  Euler  angle  transformation  from  the  inertial  to  the 
tracker  t frame  is  give  by  as  is  shomi  in  Figure  6. 
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The  llne-of -sight  (LS)  coordinate  frame  is  widely  used  in  pointing  and 

tracking  problems.  The  LS  amd  t frames  have  the  same  origin;  however, 

the  L3  frame  has  one  axis  pointing  exactly  at  the  target  while  the  t 

frame  is  misaligned  from  this  line-of -sight.  This  study  assumes  that 

the  LS  x-axis  points  directly  at  the  target  as  was  shown  in  Figure  2, 

LS 

For  perfect  tracking,  the  LS  and  t frames  are  aligned,  i.e.  = ^. 
Let  (r^g)^  denote  the  position  vector  of  the  satellite  with  respect 
to  the  tracker  expressed  in  the  inertial  frame j 


(£ts)  “ 


Transforming  this  vector  to  the  LS  frame  yields} 

(r^)“  - cf  Cr^)^  (77J 

LS  A ^ 

However,  by  definition  of  the  LS  frame  (r+  ) = 0 where  R is  the 

Lo. 

range  between  the  target  and  the  tracker.  Remembering  that 
~ G^  (for  perfect  tracking 

R cos9cos^  sin9cos^  -sin^  r^g^^ 

0 ■ -Sind  COS0  0 r^gy  (78) 

0 cosGsin^  sinGsinj^  cos0  ri.- 


-rJ^x^inG  + r^gyCosG  = 0 


so 


(80) 


Also, 


rtsx 

tan  ^ + Tf 

^tsx 


^sx  > 0 


4x<  “ 


r^gj^cosGsinJ^  + r^g^sinGsin^  + r^g^cos^  « 0 


which  implies,  after  some  algebraic  manipulation,  that 


. ^2  a/2 

^ tsx  tsy^ 

and  therefore , ,2  . 2 i /o 

fr^  + A 

/ _i  ^^tsx  ^s.V'^ 

( i2  . ii  , 'j2  a/2 

V ) * “^ts,  + nsz) 


> 0 


This  relationship  can  be  seen  in  Figure  7, 


In  practice,  perfect  tracking  in  which  the  tracker  x axis  aligns 
perfectly  with  the  LS  x-axis  will  not  be  possible.  The  misalignment 
between  the  tracker  and  LS  frames  can  be  defined  in  terms  of  two  Euler 
angle  rotations.  In  a manner  entirely  similar  to  the  previous  deriva- 
tion, the  Euler  angle  rotations  are  6v  about  the  sixis  and  6€  about 
an  intermediate  (y^  = y^^)  aucis  as  shown  in  Figure  8 and, 
r cos6vcos6€  sln6vcos6€  -sin6G'l 


-sin6v 


cosdv 


cosdvsin^ 


sin6vsin5€  cos6€ 
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It  Is  reasonable  to  assume  that  6£  and  6v  are  sufficiently  small,  even 
at  the  beginning  of  the  track,  so  the  small  angle  approximations  can 
be  used: 


W 


sin5€  ~ 5€  stn6v  ~ 6v 

cos6€  ~ 1 cos5v  ~ 1 

In  which  case,  Equation  (85)  becomes 
1 6v  -6G  ' 

-5v  1 0 

6€  0 1 


pLS 


(86) 


A solution  for  the  tine  rate  of  change  of  6€  and  6v  is  now  sought. 
Since  6€  and  6v  are  both  small  angles,  then  the  time  rate  of  change 
of  these  auigles  (6£  and  6v)  are  nothing  more  than  the  y and  z aixis 
components,  respectively,  of  the  angular  velocity  of  the  LS  frame  with 
respect  to  the  tracker  frame,  coordinatized  in  the  LS  coordinates. 
Hence, 
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1 6v  -6€  w. 


- -6V  1 0 wj 


6 0 1 I w‘ 


"iLSx  - "itx  - «ity  ^itz 


”iL3y  “ «ity  •*•  «Ux 


LS  ^ /,c  t 

.”iLSz  ■ ”ltz  " "itx 


So  since  evident  that 

” ”iLy  “ ^ity  "•  ^itx  (92) 

“ ”iLz  - '-^itz  ‘ ^itx  (93) 

^Lx  “ Vx  ^ ^Ity  ’ '^Uz  (9"*) 

Equations  (92)  and  (93)  describe  the  time  propagation  of  the  error 

misalignment  angles  6€  and  6v,  Equation  (94)  will  prove  useful  in 

the  following  development  of  the  time  evolution  of  the  line-of -sight 

LS 

emgular  velocity  vector  Wj^^S* 

In  order  to  determine  expressions  for  the  time  rate  of  change  of 
the  line-of -sight  angular  velocity  vector,  consider  the  position 
vector  of  the  satellite  with  respect  to  the  aircraft,  r^^^.  Differenti- 
ating r^g  twice  with  respect  to  time  and  applying  the  Coriolis  Theorem 
each  time  yields 
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where  the  vertical  bar  Indicates  the  frar.e  in  which  the  differentiation 
takes  place,  Coordinatizing  iDquation  (95)  in  the  Lia  frame  and 
defining  the  acceleration  of  the  satellite  relative  to  the  tracker 
along  the  line-of-slght  x,  y,  and  z axes  as 
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where  V » R = range  rate,  the  following  is  then  obtained 
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Examination  of  Equation  (99)  yields  the  following  scalar  state 


equations 
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R = V. 
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Defining  (a^a)  as  the  inertial  acceleration  of  the  satellite 
relative  to  the  tracker,  coordinatised  in  the  tracker  frame,  Equation 
(86)  can  be  used  to  obtain 
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thus,  it  follows  that 
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Equations  (lOO)  through  (102)  now  become 
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K - 4x  - 4y  - * '^£z'> 

Using  Equation  (94)  to  eliminate  from  Equations  (10?)  and  (lC8) 


yields 
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The  bracketed  £•}  terns  in  Equations  (llO)  and  (ill)  represent 
effects  due  to  the  misalignment  between  the  tracker  and  LS  frames. 

For  the  case  of  perfect  tracking  (6€  = 6v  = O)  these  terms  equal  zero. 


Measurement  2q  nations 

The  tracking  system  has  the  capability  of  measuring  the  inertial 
angular  velocity  of  the  tracker,  the  two  angular  deviations  6^  and  6v, 
and  the  range  between  the  tracker  and  the  target.  The  models  used  in 
the  measurement  equations  aire  developed  in  References  13,  14,  and  15 , 

The  inertial  angular  velocity  of  the  tracker  is  measured  in  the 
tracker  coordinate  frame  by  three  rate  gyros,  one  mounted  along  each 
tracker  axis.  While  the  system  (truth)  model  propagates  the  true  line- 
of-sight  angular  velocities  wJ^^Sy  '^^Sz  Liquations  (llO)  and  (ill)], 


measuroiaents  arc  available  only  in  the  tracker  frame  of 

and  Therefore,  raeasurcnents  of  and  considered 

by  the  filter  to  be  pseudo -racasurenents  of  and  , 

i.i^y  iLuz 

The  dominant  effects  which  contribute  to  errors  in  the  rate  gyro 
measurements  are  scale  factor  errors,  drift  errors,  g-sensitive  mass 
unbalance  errors,  misalignment  errors,  and  white  noise,  v^.  A 
suggested  gyro  rate  measurement  model  (Ref  13000)  is  given  as 
(tracker  x-axis  only) 


mx 


w^.^  + w|'.  + 2 B a.  + G„ 

itx  gsfx  itx  gmx^  i gx 


wj^.l 

-gma  -4.t 


+ V, 
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where 


mx 


"itx 


gsfx 


giax£ 


^i 


gx 


AC 

“-gma 


= measured  angular  velocity  along  x^  axis 
= true  angular  velocity  along  x^  axis 
= constant  bias  gyro  scale  factor 

= coefficients  (along  x,  y,  z directions  in  tracker  frame) 
of  the  g-sensitive  mass  unbalance  to  which  the  gyro  is 
subject 

“ accelerations  a^^^)  of  the  tracker  origin 

with  respect  to  inertial  space 
= gyro  drift  rate  along  x^  axis 

“ the  error  angle  transformation  matrix  resulting  from 
the  misalignments  of  the  three  gyros 


1*6 


AC 

“Sna 


-B 


■^gma 


21 


-3. 


B 


gma 

0 

gma. 


12  ^-"^3 


-B 


gna, 


23 


31  — 32 

i- 

“ gyro  misalignment  error  angles  between  the  desired 

gmaij 

gyro  coordinates  and  the  actual  gyro  coordinates 
(movinting  errors) 

• 3i»  i = X,  y,  z = component  of  the  vector  [ • j 
= additive  zero  mean  white  Gaussian  noise  to  account  for 
urunodeled  effects  such  as  aniso-elastic  drift,  quantiza- 
tion error,  etc. 

Suppose  a gyro  has  been  studied  in  the  laboratory.  Results  of 
the  gyro  testing  indicate  that  the  gyro  has  a drift  component  along 
the  x-axis  with  a steady  state  standard  deviation  of  a radians  per 
second  and  a process  correlation  time  of  T seconds.  These  statistics 
chciracterize  an  exponentially  time-correlated  (first  order  Harkov) 
process,  modeled  as  the  output  of  a first  order  lag  feedback  network 
driven  by  a zero-mean  white  Gaussian  noise  of  strength  Q •»  2o  /T, 

(Ref  3*4-88),  as  shown  in  Figure  9.  The  state  equation  for  0^,^  is 

% ■ "7 

The  autocorrelation  function  of  C is 

Sx 

E[Cg^(t)  Gg^(t-^T)]  = 0^ 


(113) 


The  remaining  coefficients  in  the  rate  gyro  measurement  equation 
are  modeled  as  random  biases,  (Ref  10*4-32),  as  shown  in  Figure  10. 
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With  this  model,  the  filter  is  "told"  that  the  value  of  the  variable 
does  not  change  in  time,  although  you  do  not  know  its  value  apriori. 
This  represents  a reasonable  model  for  the  remaining  coefficients 
because,  while  they  may  not  be  constant  on  a long  terra  basis,  they 
will  remain  essentially  constant  during  a ten  minute  tracking  pass. 

The  general  form  of  the  state  equation  for  these  coefficients  is 
X •=  0.  The  equation  which  describes  the  way  the  covariance  propagates 
in  time  is 


P - 0 


(114) 


This  indicates  that  the  variance  of  the  coefficient  does  not  change 
in  tine  and  that  the  initial  condition  on  P represents  the  variance 
of  the  coefficient  about  its  mean. 

The  measurements  of  the  tracker  angular  velocity  about  the  tracker 
y and  z axes  are  modeled  in  a manner  identical  to  w The  values  ixsed 

for  the  standard  deviations  in  the  process  correlation  tine  in  the 
gyro  rate  measurement  model  are  representative  of  a typical  aircraft 
rate  gyro  (Ref  13»302)» 


Quantity 

Gyro  drift 

Gyro  scale  factors 

Gyro  mass  unbalance 
coefficients 

Gyro  misalignment 
coefficients 

Additive  white 

Gaussian  noise  v^^ 


Steady  State 
Standard  Deviatipn  (g) 

-6  , 

1x10  rad/sec 

-4 

5x10 

3x10  ^ rad-sec/ra 
IxlO"^ 

1x10  ^ rad/sec 


Process 

Correlation  Time  (t) 
3600  sec 

CO 

00 

CO 

0 


j 
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The  error  misalignment  angles  6^  auid  6v  are  measured  in  the 
tracking  coordinate  frames.  Effects  which  can  degrade  these  measure- 
ments arei  deterministic  scale  factoirs,  scale  factor  errors,  angle 
track  biases,  and  angle  track  scintillation  noises.  No  attempt  has 
been  made  to  model  noises  that  are  specific  to  a typical  radar  or 
laser  ranger.  Rather,  the  measurement  model  proposed  below  is  repre- 
sentative of  a large  class  of  measurement  devices  (Ref  14: 14) 


= 6€  + S. 
m t 

(115) 

6v_  = 6v  + S^. 
m V 

+ 03p^6v  + + '’5 

(116) 

where 

6^,  6y  = true  misalignment  angles 

Sg,  =*  angle  track  scintillation  noises 
*^SFv  “ scale  factor  errors 
” angle  track  biases 

v^,  v^  = zero  mean  white  Gaussian  noises  to  account  for  unmodeled 
effects 

Both  the  angle  track  scintillation  noises  (Sg  and  3^)  and  the 
scale  factor  errors  Cgp  ) are  modeled  as  exponentially  time 

correlated  random  variables.  Scintillation  noise  is  dependent  upon 
various  factors  such  as  atmospheric  propagation,  and  amplifier 
chciracteristics  which  change  as  a function  of  time  during  a typical 
tracking  pass.  Scale  factor  errors  are  a function  of  tracker  variables 
that  undergo  a change  with  respect  to  time  (Ref  I4il5).  These  vari- 
ables are  reasonably  well  modeled  as  exponentially  time-correlated  and 
the  following  equations  are  obtained: 
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h ’ % (117) 

^ ”5  (11®) 

^SFg  “ ■^lO^SFg  '^10  (119) 

^SF^  “ *^ll°3Fv  ■*■  ”11  (120) 


where  wji^,,  w^,  v^q,  and  axe  zero  neaui  white  Gaussian  noises  with 
Qj^  = 2a^/T^,  The  represent  the  inverse  of  the  process  correlation 
time,  The  angle  track  bias  coefficients  ®ATg  and  arc 

modeled  as  random  biases  - initial  value  unloiovm  but  describable  as  a 
Gaussian  random  variable  with  mean  zero  and  a known  variance.  The 
values  used  for  the  process  standard  deviations  and  correlation  tines 
are  given  below  (Ref  4:149). 

Steady  State  Process 


Standard  Deviation  (a)  Correlation  Time 

Angle  track 
scintillations  (Sg,S^) 

1x10“^  rad 

10  sec 

Angle  measurement  , 

scale  factor  errors 
(CsFg.CsF^) 

10  " 

300  sec 

Angle  track  bias 

2x10"^  rad 

00 

(®ATg'®AT^^ 

, -6  , 

Additive  white  noise 

1x10  reid 

0 

The  model  of  the  measurement  of  reinge  is  very  similar  to  those 
of  the  angular  deviations.  Uncertainties  in  the  measurement  of  range 
are  due  primarily  to  scintillation  noise  and  bias  errors. 

*m  - “ ♦ =R  \ ''6 

where 
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range  scintillation  noise 


R = true  range 
3j^  = range  bias 

Vg  = zero  mean  white  Gaussian  noise  to  account  for  unmodeled 
effects 

The  range  scintillation  noise  is  due  to  atmospheric  effects  and 
errors  in  the  digitization  of  the  returned  signal.  The  atmospheric 


effects  in  particular  are  a direct  function  of  the  elevation  angle 
of  the  tracker  - less  scintillation  error  when  the  satellite  is 


directly  "overhead"  and  greater  errors  at  the  horizon.  The  scintil- 
lation error  will  show  a degree  of  time  correlation  during  a tracking 
pass  and  an  exponentially  time  correlated  random  variable  is  used  to 


model  this  state:  i 

Sr  = -333R  + w^  (122)  I 

where  wg  is  a zero  mean  white  Gaussian  noise  with  Q = 2a  /T  and  ^ 

33=  I/T3.  The  range  bias  is  modeled  as  a random  bias,  (initial  \ 


value  unknown,  but  describable  as  a zero  raean  Gaussian  random  variable  ^ 

j 

and  known  variance.)  Values  used  for  the  standsxd  deviations  and  j 


correlation  +jmes  are  show^.  below  (Ref  4:150). 


Quantity 

Steady  State 
Standard  Deviation  (a) 

Process 

Correlation  Time  (t) 

.i 

■i 

\ 

Range  Scintillation 

20m 

10  sec 

Range  bias 

5m 

00 

1 

Additive  white  noise 

5m 

0 

• 

1 
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1 

I 

1 

\ 

Summary  of  State  and  Measiirement  Equations 

After  augmenting  the  satellite  and  tracker  state  equations  with 
the  noise  states  needed  to  define  the  measurements,  the  truth  model 
contains  a total  of  42  states.  They  are  repeated  here  for  clarity. 


State  Equations 

(1) 

(2)  xg  - 

(3)  x^  - x^ 


Satellite  inertial  position 


W *4- V*  “dx*"l 

‘5'  ^5' V V* 

(6)  x<  “ a„„  + a„+a  +a.  +w- 

' ' o gz  mz  sz  dz  3 

ax  2 V w^^ 

/..X  . A .LS  Hsz  ''r'^lLS^ 

"LSy  “ K 


Satellite 

inertial 

velocity 


(123) 

(124) 

(125) 

(126) 

(127) 

(128) 


r"iL3y  . LS  t 
K ^ILSz  ”itx 


Hsx  ^ LS 

+ w„ 


2 v_w: 


[6v  „ - 6e  .»  1} 


/„x  • A .LS  “tsy  ''r"iLSz  LS  t 

(8)  Wl32  = WiLSz B ■ R 


R R "ILSy  "itx 

,t 

tsx  LS  Ti.  t t -1, 

- ^'iLSy  '^ity  - ’'itz^> 


(129) 

Line- 
of- 
’ sight 
angular 
velocity 

(130) 


(9)  6€  - w5:£sy  - wjty  + 6v  wf 


(10)  oC  ■ .Jfs,  - - 6e 


(11)  R - V, 


Range 


Error  misalignment 
angle 


(131) 


(132) 

(133) 


’r  - <sx  * 4.J,  - ^ 


(13)  Xj^  “ 0 Satellite  ballistic  coefficient 

(14)  Sg  - -(ii  Sg  + K4 

(15)  - ^2  Sv  + "5 


(135) 

(136) 

(137) 
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Range  scintillation 


(16) 

^ - -33  Sjj  + wg 

(17) 

V ■ A V * "7 

(18) 

°ey  ‘ "^5  °gy  * "6 

(19) 

C ■ C + w_ 

gz  6 gz  9 

(20) 

(21) 

^SFv  " “^8  ^SFy  ' 

(22) 

“0  \ 

gmxi  1 

Gyro  drift 


'11 


Angle  measurement 
scale  factors 


Coefficients  of  gyro  mass 
unbalance  (nine  eq^uations) 


(138) 

(139) 
(1^) 
(m) 

(142) 

(143) 

(1^) 


(30) 

fl  - 0 

.«“®3  J 

(145) 

(31) 

B - 0 

gmai2 

(1^) 

(32) 

B *0 

^13 

(147) 

(33) 

0 

1 

Oi 

1 

Gyro  misalignment 
coefficients  (six 

(148) 

(34) 

3 = 0 

equations) 

(149) 

(35) 

B - 0 

®^31 

(150) 

(36) 

B - 0 

(151) 

(37) 

Br-O 

Range  bias 

(152) 

(38) 

®AT6  “ ° 

Angle  track  bias 

(153) 

(39) 

BaTv  - 0 

(154) 

(40) 

^gsfx  “ ° 

1 

(155) 

(41) 

®gsfy  “ ° 

1 Gyro  scale  factors 

(156) 

(42) 

Bgsfz  “ ® 
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(157) 

Measurement 


.A 


V,  A 3'ih-Optlmal  Filter  ilodel 


Introduction 


This  chapter  presents  the  development  of  a reduced  order  system 
model  which  is  used  as  the  suboptinal  filter  model  for  the  general- 
ized ‘ionte  Carlo  Analysis  Program,  The  reduced  order  system  model 
should  be  computationally  simpler  than  the  full  system  (truth)  model. 
How  this  simplification  takes  place  is  dictated  by  the  prior  experi- 
ence, skill,  and  decisions  of  the  designer.  In  many  cases,  simplifi- 
cation is  accomplished  by  deleting  states  from  the  truth  model  that 
represent  non-dominant  effects  in  the  problem  under  study.  Typically 
this  is  accomplished  with  a corresponding  increase  in  noise  strengths 
driving  the  model  to  compensate  for  the  deleted  terms  or  states.  As 
an  example,  if  time  correlated  gyro  drift  rates  are  replaced  by  an 
additive  white  noise  in  a navigation  filter,  the  filter  performance 
nay  be  significantly  degraded.  However,  the  deletion  of  these  states 
may  not  be  noticed  for  a tracking  problem  of  ten  minutes  duration. 

In  addition  to  deleting  states,  the  designer  might  also  choose  to 
delete  terms  in  the  state  equation  that  are  about  an  order  of  magni- 
tude less  than  the  size  of  the  other  terras.  So  in  addition  to  deleting 
non-dominant  states,  non-dominant  terms  are  also  deleted.  However, 
careful  judgement  must  be  exercised  when  non-dominant  terms  are  deleted 
because  these  small  cross  coupling  terms  may  be  extremely  important  in 
meeting  system  performance  criteria. 
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The  reduced  order  filter  (a  six-state  filter)  was  suggested  by 
the  Air  Force  Avionics  Laboi-atory  to  llann  (hef  5).  The  main  idea 
behind  the  filter  is  to  use  what  is  known  about  the  dominant  forces 
acting  on  the  satellite  being  tracked  and  its  geometry  to  aid  in  the 
simplification  of  the  truth  model. 

Six  State  Filter  Equation  Development 

The  underlying  concept  for  the  development  of  the  filter  model 
proposed  in  this  section  is  to  delete  the  satellite  inertial  position 
and  velocity  states  [[x(l)  through  x(6)].  Other  information  already 
available  in  the  remaining  six  states  and  INS  data  will  then  be  used 
to  determine  the  acceleration  of  the  satellite  relative  to  the 
tracker  based  upon  the  i'u;owledge  that  the  dominant  acceleration 

of  the  satellite  can  be  described  by  a two  body  point  mass  gravity 
model.  This  approach  to  the  satellite  tracking  problem  was  suggested 
by  Air  Force  Avionics  Laboratory  personnel. 

Figure  11  represents  a typical  tracker  line-of-sight  to  satellite 
geometry.  Where 

(x^,  y^,  z^)  is  the  geocentric  equatorial  inertial  frame 
(x^^,  is  the  tracker  line -of -sight  fi'ame,  with 

pointing  at  the  satellite 

is  the  vector  from  the  center  of  the  earth  to  the  tracker 
r Is  the  vector  from  the  center  of  the  earth  to  the  satellite 

-03 

r^g  is  the  vector  from  the  origin  of  the  tracker  system  to 
the  satellite  along  the  line-of-slght  x-axis 
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Figure  11.  Tracker  Line-of-Slght  to  Satellite  Geometry 


.1 


When  the  satellite  Inertial  position  and  velocity  states  are 
deleted  from  the  system,  the  following  information  remains  available 
to  formulate  a filter  model: 

9,j^  precise  resolver  measurements  of  the  tracker  azimuth 

and  elevation  angles 
LS 

the  transformation  matrix  from  the  inertial  to  the 
llne-of-sight  frames  as  a function  of  0 and  0 
high  precision  measurements  of  the  tracker  accele- 
ration with  respect  to  inertial  space  in  tracker 
coordinates,  determined  from  the  outputs  of  three 
accelerometers  - one  mounted  along  each  of  the  tracker 
axes. 

R . w^  , , 6v  noise  corrupted  measurements  of  range, 

ra*  my'  mz*  m’  m 

tracker  angular  rates,  and  angular 
deviations  as  discussed  in  Chapter  IV. 

Consider  now  the  tracker  line-of -sight  angular  velocity  state 


equations . 


.L3  tsz 

"iLSy  = - -Tf-  - 


2 V 

^tsz  ^ ^r”iL;3y  LS 


'^iLSz  "itx 


^tsx  LS  - +. 

■r + '^iLSz  Wity  - «itzJ> 


•LS  ®tsy  ^ '^r'^LSz  LS  t 
’^iLSz  “ T H "iLSz  ”itx 

. ^-tsx  T.q  r.  4- 


+ { ^ w5;^y  [6v  wf^y  - 6€  wftz]} 


(164) 


(165) 


The  hracketted  terms  £•}  result  from  the  fact  that  the  tracker  and 


line-of -sight  frames  are  not  coincident  and  differ  by  the  two  small 

Euler  angles  6€  and  6v,  For  high  accuracy  tracking,  the  small  angular 

deviations  6€  and  6v  will  have  magnitudes  on  the  order  of  lO”'^  radians 

or  less  (Ref  10:6?).  For  the  particular  tracking  profile  used  in  this 

-11  / ? 

study,  the  bracketted  terms  are  on  the  order  of  10  radians/sec'^  while 
the  smallest  values  of  w^LSy  ^LSz  ~ radians/sec^.  Thus, 

in  the  reduced  order  filter  model  the  bracketted  terms  are  neglected 
and  replaced  by  a zero  mean  white  Gaussian  driving  noise  to  account  for 
the  increased  uncertainty  in  the  dynamics  model: 


.LS 

^•tsz 

2 V Wi'Tc„ 

X lii^y 

”iL3y 

ss  « ■ ■ 
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•L3 

_t 

^■tsy 

'Jl 

> 
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^iLSz 

“ K 

H 

^ LS  t ^ 

+ "iL3z  "itx  «1 


■^iLSy  ”itx 


+ Wr 


(166) 


(167) 


before  the  filter  was  tuned  to  give  the  best  performance,  w^  and  W2 
were  assigned  one-sigma  values  of  10  radians/sec^,  based  upon  the 
values  of  the  brac!:etted  terns  that  were  dropped. 

The  remaining  state  equations  are: 

“ ’'^y  ’ '"ity  ’4tx 

T.S  _ t 


6v  =* 

R = V 


"iLSz  ■ "itz  " 6^ 


■c2 


(169) 

(170) 
(171) 


“ Hsx  4sy  " Hsz  ^(«iLy  + ^iLSz) 

Equation  (l7l)  nay  be  simplified  using  the  same  criteria  as  for 

Equations  (l64,  I65).  For  high  accuracy  tracking,  the  terns  containing 
6€  and  6v  will  be  approximately  five  orders  of  magnitude  smaller  than 
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V^,  Therefore,  in  the  proposed  filter  model,  these  two  terms  are 
dropped  and  replaced  by  a zero  mean  white  Gaussian  driving  noise 
(w^)  of  strength  equal  to  (l  x 10“^^  m/seo^)l 


Vr  = 


tsx 


^^’^'iLSy 


w. 


LS 

iLSz 


) + 


(172) 


Consider  the  following  develcpment  of  the  acceleration  of  the 

LS 

satellite  with  respect  to  the  tracker  (a^g)  expressed  in  the  line- 
of-sip:ht  frame, 

- (St)"  + (173) 

where  the  superscript  "i"  again  indicates  coordinatization  in  the 
inertial  frame.  From  the  definition  of  the  line-of-sight  (LS)  frame 
(the  line-of-sight  x axis  - x'*^  - points  exactly  at  the  target),  we 

know  that  the  position  vector  of  the  satellite  relative  to  the  tracker 

/ vbS  LS 

(r^gj  lies  along  x . 


LS  A 


R 

0 

0 


(174) 


where  R is  the  range  between  the  tracker  origin  and  the  satellite. 


It  follows  that 
(£+=) 


LS^  i 

C.  (£+..;) 


(175) 


±^s'  -i  ''-ts-' 

From  Equation  (173)  it  is  seen  that  the  acceleration  of  the  satellite 
with  respect  to  the  tracker  expressed  in  inertial  coordinates 


(f^g)^  = (4)^  - (r^)^ 


(176) 


Elementary  astrodynamics  tells  us  that  vre  can  model  the  inertial 


acceleration  of  the  satellite  as 


r 


nie(lG)' 


(177) 


where  |jl  is  the  earth's  gravitational  constant  and  Equation  (176) 


becomes 

...  .1  ,1 

(r+ J = T~~~T^  - (^t) 


where  from  Equations  (173)  ^■nd  (175) 


(178) 


(^)^  = (Et)^  + (x^ J 


■L3 


(179) 


It  is  evident  at  this  point  that  if  the  tracker  inertia.l  position 


(r^)  were  made  available  from  an  INS  on  the  aircraft,  then  the  accele- 


ration of  the  satellite  relative  to  the  tracker  expressed  in  LS 


L3- 


coordinates  [[which  is  (a^g)  1 can  be  found  by 


/ vLS  A / '\1 

(a^.3)  = (r^g)  = (r^  J 


^ 7^'  -n 

..  i 
- (^t) 


(180) 

(181) 


.LS,.  r(^t)^'^L(£ts)^^1  p.  xi,  , X 

^LS/  \L3 

-He  (r^)  +(rts) 


Si 


(r_) 


i 3 


-(rt) 


L3 


(183) 


Under  the  assumption  that  6€  and  6v  axe  small,  the  tracker  and 
line-of -sight  frames  are  neaxly  aligned.  It  follows  that  the  accele- 


ration vector  of  the  satellite  relative  to  the  tracker  coordinatized  in 

xLS 


the  tracker  frame  - (a^g)^  - is  well  approximated  by  (a^g)  . Also, 


the  accelera.tion  of  the  tracker  with  respect  to  inertial  space  expressed 


xLS 


in  the  LS  frame  - (r^)  - is  well  approximated  by  the  inertial 
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acceleration  of  the  tracker  in  the  tracker  frame  (r^)^  which  is  derived 
from  the  outputs  of  the  accelerometers.  Equation  (I83)  can  now  be 


approximated  as 


(^3)  = — r3  - (^t) 


(184) 


Equation  (l84)  is  readily  iraplementable  as  all  quantities  in  it  axe 


available: 

L3 

is  a knovm  function  of  the  resolver  angles  ^ and  0 
(r^)^  is  provided  by  the  IKS 

(r.  ) =0  where  R is  the  range 

[0. 

/ •• 

is  derived  from  the  output  of  the  accelerometers 
(r^)*  - (r^)‘  + 4(1*3)''"  (185) 

Therefore,  the  form  of  the  s-tate  equations  for  ^i^z' 

V_,  remain  as  given  in  Equations  (I66),  (I67),  and  (l72)  with  the  compo- 


nents of  (a+  ) 


(atJ 


®tsx 

t 


(186) 


determined  by  the  corresponding  components  of  Equation  (l84). 


Filter  Measurement  Equation  Development 

The  measurement  equations  for  the  truth  model  are  summarized  at 
the  end  of  Chapter  IV.  Each  measurement  is  considered  'to  be  the  sum 
of  the  s-tate  (which  is  a realization  of  one  stochastic  process),  and 
additive  noises  (which  are  realizations  or  samples  of  other  stochastic 


processes).  For  instance,  consldei-  thCxfiteasui-ement  of  the  in,ertial 

y 

aniTulax  velocity  of  the  tracker  along  the  tracker  y cixis 


+B„w^  + a,  + C 

>ny  ity  gsfy  ity  i gy 


+ fAC  + V, 

' gma  -xt  'y  1 


(187) 


where  w^^  = t-^itx  '^ity  represents  the  true  angular  velocity 

of  the  tracker  coordinatized  in  the  tracker  frane  (it  should  be  recalled 

that  Wjjjy  and  are  considered  to  be  pseudomeasurements  of  the  line- 

LS  BS 

of  sight  angular  velocities  and  which  are  not  directly 

measurable).  The  second  through  fifth  terms  in  Equation  (18?)  are 

stochastic  models  of  the  dominant  noise  processes  that  corrupt  a rate 

gyro  measurement.  The  last  term,  Vj^,  is  a zero  mean  white  Gaussian 

noise  sequence  added  to  account  for  errors  in  the  modeling  assumptions 

and  unmodeled  higher  order  effects,  ®gsfy»  ^gmy  ’ ^gmy  ’ ^gmy 

the  elements  of  the  matrix  AG„  „ are  all  modeled  as  random  biases.  For 

-gma 

lack  of  better  information  they  are  modeled  as  zero  meain  with  a variance 
determined  from  empirical  data.  Each  of  these  stochastic  processes  is 
then  multiplied  a known  quantity  (in  the  truth  model)  - the  resulting 
product  in  e^ch  case  bcir^4  a random  process.  The  inertial  acceleration 
of  the  tracker  origin  in  the  tracker  coordinates 


(188) 


is  assumed  to  be  a deterministic  system  parameter  in  this  study.  The 

component  of  the  gyro  drift  along  the  tracker  y axis,  Cgy,  is  modeled  as 

an  exponentially  time-correlated  random  process.  In  the  filter 
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measurement  model,  it  is  assumed  that  the  total  effect  of  all  of  the 
corruptive  effects  in  each  of  the  truth  model  measurement  equations 
can  be  replaced  by  a single  zero-mean  white  Gaussian  noise  sequence. 

The  filter  measurement  model,  for  each  state,  consists  of  the  “true 
value"  plus  an  additive  white  noise  to  account  for  modeling  uncertain- 
ties, For  the  state  we  chose  as  an  exeimple, 

t t . 

w„„  = w, . + v, 

my  ity  1 

This  approach  to  modeling  the  measurements  leads  to  the  simplest  filter 
implementation.  If  the  performance  should  prove  to  be  poor  using 
this  model,  the  variance  of  v^^  could  be  increased  to  indicate  additional 
uncertainty  in  the  assumed  measurement  model.  If  performance  remains 
poor,  this  would  be  an  indication  that  some  of  the  effects  appearing 
in  the  truth  model  measurement  equations  must  be  added  - at  the  expense 
of  a higher  dimensioned  filter  - to  the  filter  measurement  equations. 

Summary  of  Filter  State  and  Measurement  Eg uations 


The  state  and  measurement  eqx;iations  for  the  reduced  order  filter 
model  are  summarized  below.  The  development  of  the  linearized  dyncimics 
and  measurement  matrices  F and  H for  the  filter  is  given  in  Appendix  B. 


VI,  Results  and  Discussion 


Introduction 

A generalized  Monte  Carlo  Analysis  Program  (i.GAp)  has  been 
developed  and  used  in  the  analysis  of  an  extended  Kalman  filter  (for 
a users*  guide  and  program  description,  see  Appendix  a).  The  purpose 
of  this  chapter  is  to  present  the  results  of  this  feasibility  study 
using  MCAP  in  evaluating  the  six-state  extended  Kalman  filter  and 
subsequently  to  discuss  and  interpret  these  results.  Prior  to 
presenting  the  results,  the  tracking  profile  used  in  the  Monte  Carlo 
analysis  and  the  philosophy  used  to  tune  the  filter  will  be  presented, 
since  they  have  direct  bearing  on  the  performance  achieved  in  this 
study. 

At  the  initialization  of  the  tracking  profile,  the  aircraft 

(tracker)  lies  on  the  Greenwich  meridian,  at  a geocentric  latitude 

of  30°  north.  For  the  duration  of  the  200  second  tracking  pass,  the 

tracker  moves  at  a constant  velocity  of  135  meters/sec  (3OO  miles/hr) 

eastward  while  maintaining  a constant  altitude  of  9000  meters  and  the 

same  geocentric  latitude.  The  satellite  is  in  a 2,0x10^  meter  circular, 

near  polar  orbit.  The  satellite  is  a relatively  small  vehicle,  with  a 

ballistic  coefficient  of  0,015  and  a solar  pressure  coefficient  equal 

to  that  of  a vehicle  with  a projected  surface  area  towards  the  sun  of 
2 

10  meter  , Initially,  the  satellite  lies  essentially  on  the  prime 
meridian  and  is  approaching  a descending  node  (descending  towards  the 
equator),  as  shown  in  Figure  12,  Pertinent  tracking  information  for 
this  flight  profile  is  summarized  in  Table  I,  The  angles  0 and  0, 


Table  I 

Typical  Tracking  Parameters 


Parameter 

Initial  Value 
(t=0) 

Final  Value 
(t=200) 

Elevation  Angle 
of  Tracker  (-^) 

56.6° 

74.5° 

Azimuth  Angle 
of  Tracker  (0  ) 

180.0° 

221.0° 

"LSy 

-4 

-9.58x10  Rad/Sec 

-4.98x10“^  Rad/Sec 

'*L3z 

2.57x10  ^ Rad/Sec 

2.47x10  Rad/Sec 

Range  (r) 

6 

2,10x10  Heters 

6.13x10^  Meters 

Range  Rate  (V^) 

3 

-7.50x10  p./Sec 

3 

-6.98x10  n/Sec 

L 
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which  define  the  coordinate  transformation  from  the  inertial  to  the 
llne-of-slght  frame  (see  Chapter  IV),  may  also  be  used  to  describe  the 
tracker  azimuth  and  elevation  angles.  The  term  is  the  elevation 
angle  and  9 the  azimuth  angle  (0°  azimuth  being  defined  as  the  condi- 
tion in  which  the  axis  is  parallel  to  the  inertial  x-axis). 

Figures  I3  and  14  present  the  tine  history  of  ^ and  9. 

The  truth  model  time  history  of  states  Wj^^y  ^L5z  x(l) 

and  x(2)3  are  given  in  Figures  15  and  I6.  Note  that  the  maximum  tracker 
angular  velocity  is  reached  at  the  end  of  the  tracking  pass.  Since 
the  range  is  decreasing,  this  condition  is  expected  from  the  geometry 
of  this  tracking  profile.  As  the  vehicle  moves,  the  tracker 
traverses  a trajectory  confined  to  a plane  perpendicular  to  the  satel- 
lite orbital  plane.  Therefore,  it  is  expected  that  the  inertial 
angular  velocity  about  the  elevation  axjs  (wj^^y)  would  continue  to 
increase  until  the  satellite  approaches  its  zenith  with  respect  to  the 
tracker.  In  ad.dition,  the  inertial  angular  velocity  about  the  azimuth 
axis  is  expected  to  continually  increase  imtil  the  satellite 

lies  in  the  plane  of  the  tracker  path  (an  azimuth  angle  of  270°). 

Thus,  while  the  tracking  profile  used  in  this  study  did  not  represent 
worst  case  conditions,  it  does  present  a highly  nonlinear  angular  rate 
history  with  which  to  evaluate  the  estimation  accuracy  of  the  Kalman 
filter. 

The  truth  model  time  history  of  states  R[x(5)  «=  range]  and 
V^x(6)  = range  rate]  are  shown  in  Figures  1?  and  18,  Note  that  the 
range  plot  is  nearly  linear  (over  the  time  period  of  interest)  and  the 
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Figure  13.  Tracker  Elevation  Angle  Time  History 
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Figure  14.  Tracker  Azimuth  Angle  Tine  History 
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Figure  15.  x(l)  State  Time  History 
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Figure  18,  x(6)  State  Time  History 
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satellite  and  tracker  are  closest  at  the  end  of  the  tracking  pass.  On 

i 

the  other  hand,  the  plot  of  range  rate  is  decreasing  with  time  and  is  : 

non-linoar  over  the  time  interval. 

Briefly  stated,  a Kalman  filter  is  considered  to  be  well  "tuned"  • 

if  the  filter  error  variance  follows  the  "true"  system  error  variance  | 

as  closely  as  possible  without  underestimating  it.  Typically  this 
will  result  in  the  lowest  possible  true  system  error  variance.  This 
is  accomplished  by  varying  the  initial  filter  covariance  matrix  I Pq  1. 
and  the  strengths  of  the  state  and  measurement  and  h]  noises  in  the 
filter  until  the  desired  "tuning"  is  achieved.  If  the  filter  under- 
estimates the  error  variance  significantly,  then  the  filter  may  diverge 
because  the  filter  is  weighting  the  filter  internal  model  too  heavily 
and  the  measurements  too  lightly.  This  is  what  Jazwinski  has  termed 
"learning  the  wrong  state  too  well"  (Ref  3.6:301-302).  On  the  other 
hand,  if  the  filter  error  variance  weights  the  measurements  too 
heavily,  then  the  filter  will  track  the  noisy  data  and  not  exploit 
the  internal  model  sufficiently.  Divergence  may  also  occur,  in  an 
extended  Kalman  filter,  if  constant  noise  strengths  (g  and  R)  are  used 
(Ref  3*7?) • Such  divergence  characteristics  can  be  remedied  to  some 
extent  by  admitting  time  varying  (such  as  piecewise  constant)  noise 
strengths  in  the  model  upon  which  the  filter  is  based.  However,  the 
scope  of  this  study  was  confined  to  achieving  adequate  traclcing  per- 
formance over  a reasonable  time  interval  with  the  use  of  a single  set 
of  noise  strengths.  Decreased  dynamic  noise  strengths  would  allow  a 
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"tighter"  tuning  during  the  initial  period  while  increased  strengths 
would  remove  (or  at  least  postpone)  the  onset  of  divergence.  In 
eventual  implementation,  tine-var:>'ing  strengths,  possibly  set  adaptively 
since  their  "best"  evaluation  would  be  trajectory  dependent,  could 
enhance  filter  performance  (Ref  3«155). 

Results  of  Filter  Tuni  n<T 

The  initial  lionte  Carlo  analysis  and  tuning  was  performed  using  j 

the  (unrealistic)  assumption  that  at  time  t the  filter  state  values 

were  equal  to  the  system  state  values  (zero  initial  error)  with  vari-  .. 

ance  equal  to  that  used  by  Mann  (Ref  5) . Even  though  this  assumption  i 

is  unrealistic  it  is  commonly  used  for  initial  tuning  of  the  (extended)  I 

1 

f Kalman  filter.  The  initial  one  sigma  values  for  Q.  and  R were  obtained  i 

I i 

using  the  apriori  estimates  shown  in  Chapter  V.  Table  II  presents  | 

these  initial  filter  tuning  parameters. 

Figures  19  tlirough  30  present  the  results  of  the  Monte  Carlo 
analysis  obtained  from  ten  simulation  runs.  In  Figures  19  through  24 
each  plot  presents  the  square  root  of  the  average  of  the  filter  covari- 
ance (denoted  by  x on  the  plot)  and  the  calcule-ted  sample  standard 

i 

deviation  of  the  true  error  (denoted  by  + on  the  plot).  Figures  25 
through  30  present  the  calculated  mean  error  (denoted  by  + on  the  plot) 
plus  and  minus  one  standard  deviation.  The  three  lines  are  not  distinct 
on  some  plots,  because  the  large  valued  mean  error  prevents  adequate 
resolution.  It  should  be  noted  that  "automatic"  scaling  is  used  to 
prepare  each  plot.  The  scale  factors  are  selected  in  order  for  the 

plot  to  fill  the  plot  area,  hence  small  valued  changes  can  be  exagge-  ; 

rated. 


77 


— ^ ^ ^ 

OV‘0  06*0  02*0  01*0  00*0 

e.Ol*  (X)jcl  ‘(  + )y0^y3  - A3a  019  (9)X 


Figure  21.  Initial  Case  - Standard  Deviation  vs  Time 
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Figure  23,  Initial  Case  - Standard  Deviation  vs  Time 
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FlKure  24.  Initial  Case  - Standard  Deviation  vs  Tims 
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Figure  25,  Initial  Case  - Mean  Error  +/“  vs  Time 
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Figure  26.  Initial  Case  - Mean  Error  +/“  vs  Time 
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Figure  27,  Initial  Case  - Mean  Error  +/-  Sg  vs  Time 
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Figure  28.  Initial  Case  - Mean  Error  +/-  vs  Time 
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Figure  29.  Initial  Case  - Ilean  Error  +/-  vs  Time 
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Figure  30*  Initial  Case  - Mean  Error  +/-  vs  Time 
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Table  II 


r 


Initial  Filter  Tuning  Parameters 


Parameter 

Matrix 

Diaigonal  Element 
Index  Number 

Value 

Units 

££0 

1 

1.0x10"^^ 

rad/ sec^ 

2 

1.0x10“^^ 

rad/sec^ 

3 

l.OxlO"® 

rad^ 

4 

l.OxlO'® 

rad^ 

5 

l.OxlO”® 

Km^ 

6 

l.Oxlo"^® 

(Km/sec)^ 

1 

1.0x10“^^ 

rad/sec 

2 

l.OxlO"^^ 

rad/sec 

3 

3.0x10"'^ 

Km/sec 

R 

1 

5.0x10"^ 

rad/sec 

2 

5.0x10"^ 

rad/sec 

3 

2.5x10'^ 

rad 

4 

2.5x10"^ 

rad 

5 

20 

m 

Figures  19  and  20  plot  the  standard  deviation  of  the  error  in  the 
estimates  of  WLSy  ^L3z  [states  x(l)  and  x(2)].  As  shovm  in  the 
plots,  the  filter  estimated  error  before  the  measurement  update  is 
greater  than  the  true  error  until  the  time  is  greater  than  180  seconds. 
At  180  seconds,  the  true  error  begins  to  diverge.  After  the  measurement 
update,  the  filter  estimated  error  is  below  the  true  error.  It  should 
be  noted  that  after  the  Initial  transient  (about  five  seconds  on  the 
plot),  the  filter  estimates  are  essentially  steady  state. 
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The  mean  error  of  state  is  plotted  in  Ficures  25  and 

26.  Since  the  filter  was  tuned  with  filter  initial  conditions  equal 
to  the  truth  model,  the  plot  begins  at  zero  mean  initial  error.  It 
should  be  noted  that  the  magnitude  of  both  plots  are  monotonically 
increasing  over  the  time  interval  and  that  the  measurement  update  tends 
tc  degrade  the  filter  estimates.  Table  III  presents  the  minimum  and 
mcLXinum  values  of  the  mean  error  starting  from  zero  initial  conditions. 

Table  III 

State  Mean  Frror  (Zero  Initial  Conditions) 


Parameter 

Minimum 

Value 

Maximum 

Value 

Units 

’^LSy  “ 

-7.3x10-7 

3.0x10"^ 

rad/sec 

«LSz  " 5^(2) 

-3.4x10"® 

2.9x10"^ 

rad/sec 

66  “ x(3) 

-7 

-1.1x10 

8.4x10“^ 

rad 

6v  “ x(4) 

-3.7x10"^ 

5.6x10"^ 

rad 

R - x(5) 

0 

2.6x10‘^3 

Km 

Vj,  “ x(6) 

0 

5.8xl0‘*'^ 

Km/sec 

Figxires  21  and  22  plot  the  standard  deviation  of  the  error  in  the 
estimates  of  the  error  states  6v  and  66  [[states  x(3)  and  x(4)[).  Note 
that  the  true  error  standard  deviation  diverges  in  both  cases.  The  mean 
error  for  states  6v  and  66  i:  plotted  in  Fj-gures  2?  and  28.  Again  note 
that  both  plots  start  from  zero  initial  error  and  increase. 

The  standard  deviation  plots  for  states  R and  V^.  [states  x(5)  and 
x(6)]  are  shown  in  Figures  23  and  24  with  the  corresponding  mean  erix>r 
plots  in  Figures  29  and  30*  Again  note  that  the  true  error  standard 


Several  Monte  Carlo  simulations  were  performed  while  varying 
values  of  ^ and  R for  each  simulation.  No  value  was  found  that  would 
prevent  the  x(5)  state  (range)  from  diverging.  Recall  from  the 
filter  equations  (Chapter  V),  the  only  x(5)  state  tuning  parameter 
is  the  measurement  noise  R^.  The  state  equation  for  the  x(5)  state 


was  modified  as  follows j 


x(5)  “ + w^ 


(201) 


where  w^  is  a "pseudonoise"  with  strength  to  be  determined  from  the 
filter  tuning.  If  a fictitious  noise,  or  "pseudonoise",  is  added  to 
the  filter  state  equations,  then  the  filter  better  accounts  for  (cannot 
neglect)  the  incompleteness  of  the  filter  representation  of  the  true 
system  (Ref  3«9-^)..  The  addition  of  this  "pseudonoise"  yields  another 
x(5)  state  tuning  parameter. 

Twelve  Monte  Carlo  simulations  (lO  runs  each)  were  then  performed 
to  tune  the  extended  Kalman  filter.  The  initial  filter  covariance 
terms  (PF^)  were  not  changed  from  the  values  shown  in  Table  II.  The 
filter  tuning  parameters  (^  and  R)  used  for  each  of  these  12  simulations 
are  shown  in  Table  IV.  The  entries  in  the  table  are  shown  in  an  abbre- 
viated form  v.here  1,0x10  is  written  1-11,  also  entries  underlined 
are  values  changed  from  the  previous  case. 

In  order  to  increase  the  standard  deviation  error  in  the  filter 

+6 

estimate  for  state  x(5),  the  "pseudonoise"  (i^^)  was  set  to  1.0x10  and 
a Monte  Carlo  simulation  performed  (case  l).  Except  for  states  x(5) 
and  x(6)  the  plots  were  the  same  as  in  the  initial  simulation  and  are 
not  shown.  Figures  31  and  32  show  the  plots  for  the  standard  deviation 
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Table  IV 


Filter  Tuning  Parameters  for  iiach 
Monte  Carlo  Simulation 


Case 

Q 

K 

1 

2 

3 

4 

1 

2 

3 

4 

5 

Initial 

1-11 

1-11 

3-7 

5-6 

5-6 

2.5-6 

2.5-6 

2.0 

1 

1-11 

1-11 

3-7 

1+6 

5-6 

5-6 

2.5-6 

2.5-6 

20 

2 

1-11 

1-11 

1-2 

1+6 

5-6 

5-6 

2.5-6 

2.5-6 

20 

3 

irii 

1-2 

1+6 

5-6 

5-6 

2.5-6 

2.5-6 

20 

4 

5-11 

5-11 

1-2 

1+6 

5-6 

5-6 

1-8 
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for  states  x(5)  and.  x(6).  For  state  x(5)  the  filter  estimate  of  the 
standard  deviation  increased  significantly  before  the  measurement  up- 
date. Also,  the  standard  deviation  of  the  true  error  has  increased. 
For  state  x(6),  the  true  and  the  filter  standard  deviation  plots  both 
increased.  The  mean  error  plots  for  state  x(5)  is  shown  in  Figure  33. 
The  mean  error  is  approximately  the  same  as  in  the  initial  simulation, 
however,  the  one  sigma  bounds  arc  wider. 
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Case  2 was  performed  by  increasins  so  that  the  error  in  state 
x(6)  would  increase.  The  only  plot  that  changed  noticeably  was  state 
x(6)  and  is  shown  in  Figure  3^.  The  true  error  standard  deviation 
remained  unchanged,  however,  the  filter  standard  deviation  increased, 

Next,  in  order  to  increase  the  error  in  states  x(l)  and  x(2), 
and  (^2  wGre  increased  for  case  3«  expected,  the  filter  estimate 
of  the  standard  deviation  for  states  x(l)  and  x(2)  increased,  so  that 
it  is  greater  than  the  true  error  standard  deviation,  which  did  not 
change  to  any  visually  apparent  degree.  These  plots  are  shovrn  in 
Figures  35  and  36.  The  other  plots  remained  the  same  and  are  not 
showwi. 

Case  k was  performed  by  decreasing  and  so  that  the  error  in 
state  x(3)  and  x(4)  would  vary  more  before  and  after  measurement  update. 

The  standard  deviation  for  states  x(3)  and  x(4)  changed  significantly, 
with  the  true  error  standard  deviation  decreasing  by  a factor  of  l/3. 

These  plots  are  shown  in  Figures  37  and  38.  All  of  the  mean  error 

plots  changed  and  are  shown  in  Figures  39  through  44,  The  most  signif-  | 

leant  change  is  noted  on  states  x(l)  through  x(4).  For  these  states  j 

I 

the  mean  error  decreased,  some  by  as  much  as  four  orders  of  magnitude,  | 

R^  and  R^^  were  decreased  again  for  case  5 ia  order  to  reduce  the  | 

filter  estimates  of  the  standai'd  deviation  for  states  3 and  4.  - 

Again,  as  shown  in  Figures  45  and  46,  the  standard  deviation  for 
states  x(3)  and  x(4)  changed  significantly.  Also  all  of  the  mean 

error  plots  changed  noticeably  and  are  shown  in  Figures  4?  through  52.  - 
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Case  6 was  performed  by  increasinf^  so  that  the  error  in  state 
x(5)  would  ' ry  less  before  and  after  measurement  update.  All 
standard  deviation  and  mean  plots  changed  significantly  and  are  pre- 
sented in  Figures  53  through  64.  host  noticeable  is  the  decrease  in 
the  true  error  standard  deviation  on  all  plots,  some  by  as  much  as  two 
orders  of  magnitude. 

In  order  to  reduce  the  filter  standard  deviation  estimate  for 
state  x(5)  Qij,  was  decreased  for  case  7.  The  standard  deviation  for 
states  x(5)  and  x(6)  changed  and  are  presented  in  Figures  65  and  66, 
The  most  noticeable  change  in  Figure  65  is  the  decrease  in  the  vari- 
ations of  the  filter  estimates  of  the  standard  deviation  before  and 
after  the  measurement  update.  In  Figure  66  the  standard  deviation 
for  the  true  error  decreased.  Figures  6?  through  72  present  the  mean 
error  plots. 

Case  8 was  performed  by  decreasing  the  value  of  and  so  that 
the  filter  standard  deviation  estimate  for  state  x(5)  would  be  reduced 
further.  The  filter  estimate  of  the  standard  deviation  for  both 
states  x(5)  and  x(6)  decreased.  However  for  both  states,  the  true 
error  standard  deviation  increased  significantly  as  noted  in  Figures 
73  and  74.  The  mean  error  plots  are  presented  in  Figures  75  through 
80. 

Case  8 yielded  too  much  reduction  in  the  filter  standard  deviation 
estimate  for  state  x(5)  so  vras  increased  by  an  order  of  magnitude 
to  perform  case  9.  The  filter  estimate  of  the  standard  deviation  for 
state  x(5)  increased  slightly,  however,  the  true  error  standard  devi- 
ation decreased  by  a factor  of  l/2,  as  shown  in  Figure  81,  For  state 
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x(6),  the  filter  estimate  remained  the  same,  but  the  true  error 
decreased  by  a factor  of  l/2  as  shown  in  Fif^ure  82.  Figures  83  and 
84  present  the  mean  error  plots  for  these  two  states. 

Case  10  was  performed  by  decreasing  by  an  order  of  magnitude 
so  that  the  filter  standard  deviation  estimate  for  state  x(5)  would 

i 

be  reduced  further.  Both  the  filter  estimate  and  the  true  error 
standard  deviation  decreased  for  state  x(5)  and  is  shown  in  Figure 

85.  Again,  the  true  error  decreased  for  state  x(6)  as  shown  in  Figure 

86.  Figures  8?  and  88  present  the  mean  error  plots  for  these  two 
states. 

Since  the  true  standard  deviation  was  greater  than  the  filter 
estimate,  R was  increased  by  a factor  of  five  to  perform  case  11. 

5 

For  state  x(5),  the  filter  estimate  of  the  standard  deviation  increased 
and  the  true  error  decreased  as  shown  in  Figure  89.  Figure  90  shows 
that  the  true  error  also  decreased  for  state  x(6).  The  mean  error 
plot  for  state  x(5)  Is  shown  in  Figure  91. 

Case  12  was  performed  by  decreasing  R^  by  approximately  l/2  so 
that  the  filter  estimate  would  be  closer  to  the  true  standard  deviation 
at  the  end  of  the  simulation  time.  The  filter  estimate  of  the  standard 
deviation  decreased  and  the  true  error  increased.  The  complete  results 
for  this  case  are  shown  in  Figures  92  through  103,  Since  this  is  only 
a feasibility  study,  further  tuning  was  not  attempted. 

These  results  demonstrate  that  all  filter  state  estimates  of  th 
true  standard  deviation  error  are  "tunable".  As  the  filter  standard 
deviation  is  changed  for  some  of  the  states,  cross-coupling  can  affect 
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Figure  101,  Case  12  - Mean  Error  +/-  s vs  Time 
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the  estimates  for  the  other  states  significantly.  Most  noticeable 
in  these  plots  is  that  the  filter  estimate  of  the  mean  error  is 
biased.  This  bias  appears  to  grow  without  bound  for  the  time  period 
of  interest.  The  next  chapter  will  present  some  recommendations  for 
compensating  this  bias. 

Extended  Simulation  Runs 

Additional  Monte  Carlo  simulations  were  performed,  using  the 
parameters  for  case  12,  to  obtain  the  results  for  5»  20,  and  40 
simulation  runs.  These  additional  runs  were  performed  in  order  to 
show  the  change  in  calculated  statistics  with  additional  data  points 
(samples).  For  each  of  these  simulations  (case  A - 5 runs,  case  B - 
20  nms,  and  case  C - 40  runs),  the  average  filter  estimate  of  the 
standaxd  deviation  for  each  corresponding  state  was  essentially  the 
same.  However,  there  is  some  variation  in  the  plots  of  the  true  error 
standard  deviation  between  the  cases.  Most  noticeable  is  the  true 
error  for  state  x(5).  Figures  104  through  106  show  the  standard  devi- 
ation plots  for  state  x(5)s  case  A,  case  B,  and  case  C,  The  mean 
error  plots  change  only  slightly  from  case  to  case  with  the  number  of 
simulation  runs  and  are  not  shown. 

These  additional  simulation  runs  were  performed  because  Bucy, 
Hecht,  and  Senne  (Ref  17:79-85)  make  the  statement  that  significant 
conclusions  (having  a statistical  confidence)  cannot  be  made  on  the 
basis  of  only  ten  simulation  runs.  However,  these  results  indicate 
that  if  a sufficient  number  of  simulation  runs  are  performed  to  obtain 
reasonable  confidence  in  the  resultant  statistics,  that  the  results 
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will  be  similar  to  those  presented  in  Figures  92  through  103  (for 
ten  simulation  runs).  This  result  indicates  that  the  filter  can  be 
tuned  using  only  a relatively  small  number  of  simulation  runs.  Once 
the  desired  tuning  is  achieved,  the  number  of  runs  can  be  increased  to 
gain  confidence  in  the  resultant  statistics. 

Recovery  from  Initial  Condition  Error 

While  it  is  generally  accepted  practice  to  initially  tune  the 
(extended)  Kalman  filter  with  zero  initial  condition  error  and  with 
Pq  as  the  best  apriori  knowledge  of  the  initial  condition  error,  this 
is  not  the  conclusion  of  the  filter  performance  analysis.  Since,  in 
the  real  world,  the  starting  conditions  of  the  filter  will  almost 
always  be  different  from  the  "true"  system,  the  filter  should  be 
analyzed  with  non-zero  initial  error.  The  question  then  arises  as  to 
what  initial  error  should  be  used.  One  accepted  method  is  to  displace 
the  filter  initial  condition  from  the  truth  model  by  one  sigma  of  the 
aprioi'i  standard  deviation  and  observe  the  change  in  the  mean  error 
plots  for  a Kalman  filter  and,  for  the  extended  Kalman  filter,  the 
standard  deviation  plots  as  well. 

One  Monte  Carlo  simulation  (lO  runs)  was  performed  with  one 
filter  state  having  an  initial  condition  error  displaced  by  one  sigma 
of  the  apriori  standard  deviation.  Repeating  for  each  state  variable 
Initially  offset,  the  standard  deviation  plots  were  the  same  as 
Figures  92  through  97.  However,  the  mean  plots  did  change  slightly 
and  are  discussed  below. 
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Figures  10?  and  108  present  the  mean  error  plots  for  states  x(l) 
and  x(4)  when  the  simulation  was  performed  with  x(l)  having  an  Initial 
condition  error.  Note  In  Figure  10?  at  time  t^  state  x(l)  starts 
viith  an  initial  error  of  -3.2x10  ^ rad/sec  and  recovers  to  the  plot  of 
Figure  98  within  10  seconds.  This  initial  transient  is  also  noted  in 
state  x(4),  Figure  108,  v;hen  compared  to  Figure  101. 

The  Monte  Carlo  simulation  was  subsequently  performed  with  an 
x(2)  initial  condition  error- and  the  mean  error  plots  for  states 
x(2)  and  x(3)  are  presented  in  Figures  109  and  110.  In  this  case, 
state  x(2)  and  x(4)  have  an  initial  transient  when  compared  to 
Figures  99  and  100. 

Figures  111  and  112  present  the  mean  error  plots  for  states  x(3) 
and  x(2)  when  the  simulation  was  performed  with  x(3)  having  an  initial 
condition  error.  In  this  case,  x(3)  has  an  initial  transient  and 
recovers  to  Figure  100  within  10  seconds.  This  initial  transient  is 
also  noted  in  state  x(2);  Figure  111  when  compared  with  Figure  99* 
however,  the  transient  lasts  15  seconds. 

For  the  case  with  x(4)  initial  condition  error,  similar  results 
were  observer’,  with  state  x(4)  having  the  transient  and  the  transient 
also  being  observed  in  state  x(l).  When  the  ilonte  Carlo  simulation 
was  performed  with  states  x(5)  and  x(6)  having  an  initial  condition 
error,  no  change  was  noted  in  the  mean  error  plots.  These  results  are 
not  presented. 

In  summary,  if  states  x(l)  or  x(2)  have  an  initial  condition  error, 
then  this  will  result  in  an  initial  transient  in  states  x(3)  and  x(4) 


or  x(2)  and  x(3)  respectively.  This  transient  recovered  to  the  case 
of  zero  initial  condition  error  within  10  seconds.  This  result  is 
also  noted  if  the  initial  condition  is  reversed;  x(3)  or  x(4)  viith 
initial  condition  error  and  the  transient  in  states  x(2)  or  x(l).  The 
transient  for  this  second  set  of  cases  lasted  15  seconds. 

Random  Initial  Condition  Error 

In  order  to  simulate  a condition  where  the  filter  initial  condi- 
tions differ  from  the  true  system  values,  the  filter  values  assamed 
initial  values  which  were  changed  from  the  true  value  by  an  additive 
white,  Gaussian  noise  having  strength  eq,ual  to  the  initial  covariance. 
This  was  accomplished  by  calling  subroutine  HCISS  (see  Appendix  A) 
in  the  first  section  of  subroutine  XFDOT.  The  value  for  the  variable 
1 RilS  was  the  sc[uare  root  of  the  initial  covariance  diagonal  element. 

^ The  result  of  IJ0I3S,  XNOISS,  was  added  to  the  initial  filter  state 

value  yielding  a random  initial  condition  for  the  initial  filter 

■ 

states  with  every  simulation  run.  While  varying  the  truth  model 
states  would  be  more  realistic,  the  initial  condition  value  would  be 

, more  difficult  to  obtain  because  of  the  interrelationship  in  the  sys- 

1 

tern  states,  (a  variation  in  satellite  position  affects  Wj^gy,  Wj^-g* 

R,  and  V^).  Varying  all  the  filter  initial  condition  values  roughly 
corresponds  to  the  truth  model  starting  from  a different  initial  condi- 
tion for  each  run. 

The  results  of  this  iionte  Carlo  simulation  (lO  runs)  is  shown  in 
Figures  113  through  124.  The  average  filter  estimate  of  the  standard 
deviation  is  the  same  as  case  12  (Figures  92  through  9?),  however,  the 
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true  error  standard  deviation  has  changed  for  all  states  as  expected. 
Note  that  in  all  cases  the  standard  deviation  for  the  filter  estimate 
is  greater  than  the  true  error  so  the  filter  is  liot  divergent.  The 
mean  error  plots  have  also  changed  when  compared  with  Figures  98 
through  103,  The  most  noticeable  result  in  states  x(l)  through  x(4) 
is  an  initial  standard  deviation  transient  lasting  approximately  I5 
seconds. 

These  results  compare  with  the  previous  section,  where  now  Instead 
of  only  one  filter  state  having  an  initial  condition  error,  all  states 
have  an  initial  error.  The  transient  noted  is  different  (a  standard 
deviation,  no*  i mean),  however,  the  transient  duration  is  approximately 
the  same  (15  seconds). 

Comparison  with  Previous  Results 

As  noted  in  the  users'  guide  for  MGAP  (see  Appendix  A)  one  of 
the  stringent  assumptions  made  for  the  GCAI-'  covariance  analysis  is  ^ 
that  total  impulsive  control  is  applied  to  all  states.  The  problem 
investigated  in  this  study  (and  presumably  in  the  study  of  Kann  and 
Mitchell)  applied  impulsive  feedback  to  only  two  of  the  six  filter 
states.  Since  the  covariance  analysis  used  in  the  previous  studies 
(GCAP)  assumed  feedback  of  all  error  states  to  the  system  and  several 
Fortran  coding  errors  were  detected,  no  direct  comparison  can  be  made 
with  the  previous  studies.  Consequently,  the  results  of  this  study 
Invalidate  the  results  of  their  previous  work.  However,  it  should  be 
noted  that  most  of  the  corrected  computer  subroutines  developed  by 
Mann  were  easily  adaptable  to  the  requirements  of  the  generalized 
Monte  Carlo  Analysis  Program  (MCAP), 
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Additional  Comments 

The  tuning  of  a (linear  or  extended)  Kalman  filter  is  a time 
consuming,  sometimes  tedious  and  frustrating  experience.  In  particu- 
lar, if  the  initial  tuning  parameters  are  not  "close"  to  a (proper 
but  unknown)  range  of  values,  then  the  computations  made  in  the  finite 
wordlength  computer  may  exceed  the  machine  limits  ("blow  up"). 

If  this  occurs  then  the  analyst  must  modify  the  simulation  (through 
insight  and  experience)  in  order  to  determine  an  acceptable  starting 
set  of  values.  In  investigating  this  problem  three  program  modifica- 
tions were  used  to  start  the  simulation: 

1.  The  filter  and  system  state  equations  were  decoupled  (no 
feedback),  the  first  measurement  update  time  set  to  the  end 
of  the  problem,  and  the  filter  covariance  propagation  equa- 
tions bypassed.  These  modifications  allowed  the  true  error 
standard  deviation  transients  to  reach  steady  state.  Filter 

values  were  then  modified  considering  these  values. 

2.  In  this  problem  formulation,  no  measurement  is  available  to 
the  filter  for  the  range  rate  state  = x(6)j.  When  the 
filter  performed  the  measurement  update,  the  Kalman  filter 
gain  matrix  K would  cause  the  filter  estimate  of  range  rate 
to  be  grossly  in  error.  This  error  quickly  rippled  through 
the  filter  state  equations,  and  the  simulation  exceeded  the 
machine  numeric  capabilities.  The  program  modification  was 
to  set  the  filter  estimates  of  range  rate  equal  to  the  true 
value  after  the  update.  After  several  simulations,  the 
appropriate  tuning  parameters  were  determined  which  gave 

196 


r ’ — ' " ^ ^ 

better  estimates  of  the  range  rate  and  so  the  modification 
was  no  longer  needed, 

3.  The  idea  of  a transformation  of  state  variables  is  used  in 
this  program  modification.  Since  this  problem  did  not  use 
error  states  for  the  filter  or  system  equations,  calcula- 
tions are  performed  using  large  and  small  state  values, 
causing  severe  numerical  problems  on  a finite  wordiength 
computer.  If  the  large  state  values  are  scaled  down,  then 
this  problem  is  eliminated.  To  better  understand  this  idea, 
the  development  is  presented  in  the  next  paragraph. 

The  idea  of  rescaling  the  problem  is  best  understood  by  considex'ing 
the  linear  filter  state  and  measurement  equations: 

X = F X + G w (202) 

z = H X + V (203} 

How,  it  is  desired  to  rescale  some  of  the  states  (the  large  valued 
states)  to  have  values  closer  to  the  other  states.  If  the  rescaled 
state  vector  is  denoted  by  x*,  then  the  following  relationships  are 
fo\md: 


x'  = T X 

(204) 

F*  = T F t”^ 

(205) 

G ’ = T G 

(206) 

H • = H T“^ 

(20?) 

T 

P’  = T P T 

(208) 

where  T is  the  diagonal  transformation  (rescaling)  matrix  that  con- 
tains the  scale  factors  along  the  diagonal.  In  the  case  of  ncn-linear 
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state  and  measurement  equations,  then  the  relationships  are  best 
found  by  writing  out  the  matrix  equations  to  determine  the  relation- 
ships. 

The  problem  investigated  in  this  study  used  this  idea  of  matrix 
transformation  to  rescale  the  filter  range  and  range  rate  states 
[x(5)  and  x(6)]  and  the  satellite  position,  satellite  velocity,  range 
and  range  rate  states  [x(l)  through  x(6),  x(ll),  and  x(l2)]  in  the 
truth  model. 
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V 

VII.  Conclusions  and  Recommendations  f 

~ ' I 

Conclusions 

A generalized  iionte  Carlo  Analysis  trogram  (HCAP)  which  evalu-  | 

ii 

ates  the  performance  of  either  the  linear  or  extended  Kalman  filter 
has  been  developed.  MCAP  Is  similar  in  structure  to  a widely  used 
General  Covariance  Analysis  Program  (GCAP)  LRef  2],  and  so  users  should 
have  little  difficulty  performing  a Monte  Carlo  analysis  after  the 
(linearized)  Kalman  filter  has  been  analyzed  by  a covariance  analysis. 

Also,  with  MCAP,  the  user  can  now  analyze  an  extended  Kalman  filter 
properly  because  the  state  values  are  available  in  the  Monte  Carlo 
analysis.  As  noted  in  Chapter  VI,  MCAP  can  be  used  to  evaluate  the 
(extended)  Kalman  filter  under  a variety  of  conditions,  for  example, 
results  of  performance;  to  zero  initial  condition  error,  to  single/ 
multiple  state  initial  condition  error,  or  to  random  initial  condition 
error.  Thus,  MCAP  is  a versatile  and  viseful  computer  program  for 
the  design  and  analysis  of  the  (extended)  KaLman  filter. 

The  feasibility  of  using  a six  state  extended  Kalman  filter  for 

!‘ 

estimating  the  states  of  the  tracking  system  (^LSy*  i 

and  Vj.)  has  oeen  demonstrated.  As  noted  in  Figures  19  through  103,  j 

when  the  tuning  parameters  are  varied,  the  filter  performance  changes 
accordingly,  and  so  the  filter  can  be  tuned  for  optimum  performance. 

In  general  one  can  postpone  the  onset  of  divergence  with  constant  ^ 

parameters,  however  better  performance  can  be  achieved  with  nonstation-  v 

i 

ary  noise  or  using  an  adaptive  filter.  It  is  also  noted  that  the 
filter  produces  biased  estimates.  These  biased  estimates  are  very 
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characteristic  of  extended  Kalman  filters,  and  is  due  to  the 
neglected  higher  order  effects  inherent  in  using  linear  perturbation 
techniques.  The  more  pronounced  the  nonlinearities  are  in  a given 
application,  the  more  serious  one  can  expect  performance  to  be 
degraded  by  this  effect.  In  order  to  compensate  for  these  bias  errors, 
bias  correction  terms  can  be  added  to  the  filter  model  (Ref 

As  noted  in  Chapter  III,  a filter  often  is  considered  to  be 
tuned  when  the  filter  estimated  standard  deviation  deviations  are 
approximately  equal  to  the  computed  standard  deviations  or  when  the 
true  error  standard  deviations  are  smallest  in  some  sense.  This 
"covariance  matching"  technique  may  not  alx-fays  yield  the  best  perform- 
ance for  a given  set  of  conditions.  This  is  observed  for  the  case  of 
zero  initial  condition  error  by  comparing  case  7 (Figures  6?  through 
72)  and  case  12  (Figures  93  through  IO3)  mean  error  plots j Sven 
though  case  12  is  "more  appropriately"  tuned  (compare  states  5 and  6) 
better  performance  is  indicated  in  case  7 by  the  mean  error  plots. 
Consequently,  the  filter  should  be  tuned  1:^  considering  the  problem 
conditions  (how  accurately  the  initial  conditions  are  known,  typical 
versus  worst  case  tracking  profiles,  etc,)  and  obtaining  the  best 
performance  in  both  the  standard  deviation  plots  and  the  mean  error 
plots. 

Recommendations 

Even  though  IICAP,  as  written,  is  a general  Monte  Carlo  analysis 
program,  further  improvements  can  be  added.  These  improvements  will 
be  discovered  after  several  different  estimation  problems  are 
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analyzed  using  liCAP.  Two  additions  currently  come  to  mind.  The 
first  addition  concerns  the  white  Gaussian  noise  generator  presently 
Implemented  in  subroutine  NCI3S.  This  noise  generator  executes 
rapidly  and  provides  random  numbers  which  are  approximately  white 
and  Gaussian.  If  a better  but  slower  executing  noise  generator 
is  desired}  the  implementation  shown  in  Ref  1?  can  be  substituted  for 
NOISE.  In  addition  to  providing  better  white  Gaussian  numbers,  the 
implementation  shown  in  Ref  1?  features  computer  to  computer 
(different  types  of  computers)  repeatibility. 

Since,  in  general,  the  user  may  want  to  specify  the  filter  method 
of  integration,  the  second  improvement  would  be  to  modify  HCAP  to 
add  this  capability.  Several  filter  integration  methods  could  be 
programmed  and  the  user  would  select  the  desired  method  with  a problern 
dependent  control  parameter  ($U3RCTL), 

Further  woi'k  can  also  be  recommended  for  the  filter  analyzed  in 
this  study.  The  performance  variation  of  the  filter  should  be 
studied  for  many  different  tracking  profiles.  The  filter  can  then  be 
tuned  (both  standard  deviation  and  mean)  consideriiig  these  several 
trajectories.  This  additional  information  will  give  insight  into 
the  nature  of  the  bias  correction  terms  needed  to  improve  filter 
performance.  These  bias  correction  terms  should  then  be  added  to 
the  filter  model.  If  these  correction  terms  still  do  not  improve 
performance  adequately,  then  either:  l)  a time  varying  and  R as 
some  unknown  function  of  the  angular  rates  ajid/or  range  data  could 
be  determined  after  studying  several  representative  trajectories. 
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or  2)  an  adaptive  ^ technique  could  be  used  based  upon  the  real  time 
evaluation  of  the  innovations  (residuals)  sequence(s). 

Several  simplifying  assumptions  were  made  in  order  to  limit  the 
scope  of  this  study.  Additional  analysis  can  be  performed  by 
relaxing  these  assumptions.  In  particular:  The  resolver  measure- 
ments of  0 and  0 were  considered  to  be  perfect.  These  paurameters 
should  be  modeled  stochastically  and  included  in  a future  filter 
analysis  to  determine  the  performance  sensitivity  of  the  filter  to 
these  variables.  Also,  the  performance  sensitivity  of  the  filter 
to  accelerometer  noise  corruptions  should  be  studied. 
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USERS'  GUIDE  FORj 

A GENERALIZED  MONTE  CARLO  Al^ALYSIS  PROGRAM  FOR  KALi’iAN  FILTER  DESIGN 


Section  I 
Introduction 

Overview 

A large  class  of  estimation  problems  is  concerned  with  finding  an 
optimal  estimate  of  some  quantity  (an  unknown  parameter,  a random 
variable,  or  a random  signal)  given  noise-corrupted  measurements  of 
a function  of  this  quantity  (corrupted  by  noise).  In  particular,  there 
are  many  problems  with  aerospace  applications,  such  as  navigation, 
guidance,  and  weapon  delivery  systems,  which  require  that  noise- 
corrupted  measurements  be  used  to  estimate  certain  physical  parameters 
precisely.  For  instance,  the  techniques  used  to  combine  these  external 
measurements  with  Inertial  Navigation  System  (iKS)  outputs  fall  into 
two  general  categories:  conventional  continuous-feedback  damping  and 
Kalman  filter  damping  (Ref  l;l).  The  trend  in  recent  years  has  been 
toward  extensive  use  of  Kalman  filter  techniques. 

In  complex  aerospace  systems,  the  number  of  parameters  needed  to 
describe  the  system  (the  true  system  model)  accurately  may  be  extremely 
large.  The  real  time  implementation  of  this  true  system  model  within 
the  "optimal"  Kalman  filter  is  often  not  practical  because  of  the 
large  number  of  parameters,  memory  requirements,  and  resultant  compu- 
tational burden  Imposed  on  the  airborne  computer.  In  order  to  obtain  a 
Kalman  filter  which  is  computationally  feasible,  intentional  modeling 
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approximations  are  Introduced  by  deleting  nondominant  states  and 
terms  from  the  system  model.  Thus  a "suboptlmal"  Kalman  filter  based 
on  this  reduced  order  model  will  give  poorer  perfcnnance,  from  a 
statistical  standpoint,  than  the  optimal  filter.  Cc.'puter  simula- 
tion techniques  are  used  to  develop  the  required  compensation  (filter 
tuning)  and  to  study  the  effects  of  uncertainties  (an  error  analysis) 
in  the  true  system  model.  Two  types  of  simulation  techniques  are 
commonly  used:  the  covariance  analysis  and  the  Monte  Carlo  analysis. 

A covariance  analysis  generates  the  second  order  statistics  of 
the  error  between  suboptlmal  Kalman  filter  state  estimates  and  the 
corresponding  parameters  from  the  true  system  model.  A General 
Covariance  Analysis  Program  (GCAP)  has  been  developed  by  the  United 
States  Air  Force  Avionics  Laboratory  and  is  coming  into  widespread 
use  (Ref  2),  The  covariance  analysis  is  somewhat  limited  because 
stringent  assumptions  (see  Appendix  a)  are  necessary  for  the  analysis 
results  to  be  a valid  depiction  of  the  error  characteristics. 

The  Monte  Carlo  analysis,  on  the  other  hand,  actually  conducts  a 
sample-by-sample  simulation  using  random  number  generators  and 
shaping  filte-^s  to  generate  the  random  error  sources.  The  noisy  system 
measurements  are  then  processed  by  the  suboptlmal  (extended)  Kalman 
filter  algorithm  to  generate  the  filter  estimates.  If  many  of  these 
simulations  are  conducted,  then  the  statistics  of  the  error  between 
the  filter  estimate  and  the  truth  model  can  be  computed.  Because  of 
the  larger  number  of  simulations  necessary  to  determine  the  performance 
of  one  set  of  Kalman  filter  parameters,  a ilonte  Carlo  analysis  can 
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Consequently  when  a covariance  analysis  can  be  performed,  the  Honte 
Carlo  analysis  Is  viewed  as  a step  to  be  performed  after  a covariance 
analysis.  The  iionte  Carlo  analysis  is  thus  used  to  fine  tune  the 
Kalman  _*lter  and  validate  the  filter  performance  after  the  covariance 
analysis  has  been  performed. 

Since  a Monte  Caxlo  analysis  of  promising  filter  designs  should 
be  performed  after  the  covariance  analysis,  there  is  a definite  need 
for  a general  iionte  Carlo  Analysis  Program  (iiCAP).  This  program 
should  be  similar  in  structure  to  GGAP  so  that  users  can  easily 
transition  from  a covariance  analysis  to  a Monte  Carlo  analysis  of  the 
Kalman  filter  under  study.  It  is  also  desirable  that  the  problem 
dependent  computer  subroutines  be  applicable  to  both  GCAP  and  liCAP, 

This  report  presents  a description  and  users  guide  for  the  general- 
ized Monte  Carlo  Analysis  Program  (iiCAp).  Section  II  presents  the 
Kalman  filter  and  extended  Kalman  filter  equations  used  in  the  Monte 
Carlo  analysis  of  the  suhoptimal  Kalman  filter.  Section  III  briefly 
describes  the  Monte  Carlo  analysis  simulation  method  and  how  the  error 
statistics  are  computed.  Section  IV  describes  the  general  Monte  Carlo 
Analysis  Program  (MCAP)  structure  and  provides  a brief  description  of 
the  constituent  subroutines.  The  documentation  presented  is  sufficient 
to  allow  the  reader  to  use  MCAP  for  the  analysis,  evaluation,  and 
design  of  suboptimal  Kalman  filters.  This  report  is  concluded  with  a 
sample  problem  using  MCA?,  presented  in  Section  V, 
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Section  II 

Extended  Kalman  Filter  Equations 

Introduction 

This  section  presents  the  propagation  and  update  equations  for 
both  linear  and  extended  Kalman  filters.  A Kalman  filter  Is,  concisely 
stated,  an  optimal  recursive  data  processing  algorithm  for  the  deter- 
mination of  the  states  or  parameters  of  a system  using  noise  corrupted 
measurements  (Ref  6tl6).  If  a physical  system  of  interest  can  be 
modeled  by  a set  of  ordinary  linear  differential  equations  emd  linear 
measurements  with  system  auid  measurement  noises  which  are  white  and 
Gaussian,  then  the  Kalman  filter  will  provide  the  best  estimate  of 
the  system  states.  However,  In  many  cases  of  practical  interest, 
physical  systems  must  be  represented  by  a nonlinear  set  of  differential 
equations  and/or  nonlinear  measiirement  equations.  For  such  problems, 
it  is  often  convenient  to  linearize  the  system  equations  about  some 
assumed  set  of  nominal  conditions  and  use  developed  algorithms  (such 
as  the  extended  Kalman  filter  which  uses  reevaluation  of  the  nominal 
at  each  mea5\u:ement  time)  for  estimation  about  these  nominal  conditions 
(Ref  7 I 57). 

Considering  the  approximations  necessary  and  the  fact  that  there 
is  no  "best"  suboptimal  filter,  the  extended  Kalman  filter  gain  and 
covariance  propagation  equations  have  the  same  form  as  the  Kalman 
j filter  equations,  but  are  linearized  about  the  current  state  estimates. 

The  linearization  is  a first  term  approximation  to  a Taylor  series  ij 

ij 

^ exi)ansion  .to  the  state  estimate.  Higher  order  and  more  exact  approxi-  |j 

I 

mations  can  be  achieved  by  using  more  terms  of  the  Taylor  series  | 

expansion  for  the  nonlinearities  and  by  deriving  approximate  recursive  f 
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relations  for  the  higher  moments  of  the  state  vector  x (Ref  8j184), 
yielding  higher  order  nonlinear  filters.  As  might  be  expected,  if 
the  system  nonlinearities  are  significant,  then  neglecting  the  higher 
order  terms  will  result  in  biased  estimates.  However,  when  compared 
to  the  extended  Kalman  filter,  the  higher  order  filters  are  both  more 
complex  and  more  costly  in  terms  of  computer  implementation.  For  this 
reason,  the  extended  Kalman  filter  is  often  considered  first  in  non- 
linear estimation  problems. 

Notation 

This  study  has  adopted  the  notation  presented  by  Wrigley,  Hollister, 
and  Denhard  (Ref  9*20-23).  A vector  (represented  by  a letter  with  an 
underbar,  x)  is  considered  to  be  a geometric  entity  in  real,  three 
dimensional  space.  The  vector  represents  some  physical  quantity  which 
has  both  magnitude  and  direction.  When  the  physical  quantity  is 
measured  with  respect  to  a coordinate  frame,  the  vector  is  said  to  be 
coordinatizod  in  that  frame.  The  three  nismbers  associated  with  the 

I 

i mathematical  vector  axe  the  components  of  the  physical  vector  relative 

I 

j to  the  specified  coordinate  frame.  As  an  example,  if  x is  coordinatized 

! . 

in  the  "i"  frame,  the  vector  would  be  denoted  by  x^,  a three- tuple  of 

numbers.  Another  vector  of  interest  is  the  physical  angular-velocity 

vector,  generally  denoted  by  a w.  The  angular  velocity  vector  will 

have  two  subscripts,  as  an  example  w^^.  The  subscripts  indicate  that 

the  angular  velocity,  w,  is  cf  frame  b with  respect  to  frame  n, 

coordinatized  in  the  1 frame. 
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A 3*3  direction  cosine  matrix  is  used  in  this  study  to  transform 
the  components  of  a vector  in  one  frame  to  those  in  another.  The 
direction  cosine  matrix  is  represented  hy  a capital  and  an  underbar. 
Associated  with  the  letter  C will  be  a subscript  for  the  frame  from 
which  the  transformation  is  made,  and  a superscript  for  the  new  coordi- 
nate frame.  As  aji  example,  would  transform  the  vector  x from  the 
b frame  to  the  n frame,  as  x*^  = x^. 

Where  it  is  necessary  to  address  individual  components  of  a vector 
coordinatized  in  a specific  frame,  the  vector  will  be  specified  and 
subscripts  used  to  indicate  individual  components.  For  example; 


where  T denotes  the  transpose  operator. 


Definitions 

Listed  below  are  some  of  the  definitions  used  in  this  chapter; 
xCt^)  = system  state  at  time  tj^  (n-vector) 

x(tp  = filter  estimate  prior  to  incorporating  a measurement  at 
time  tj^  (n-vector) 

x(t^)  = filter  estimate  after  incorporating  a measurement  at  time 
tj^  (n-vector) 

^(tj^)  = m vector  of  measurements  at  time  tj^ 

w(tj^)  = system  dynamics  white  Gaussiaji  noise  s-vector,  independent 

of  x(tjj),  where  5Cw(t)]  = 0,  and  i£w(t)w(T)'n  = S( 't)^  ( t-T) , 

[[w(t)  is  assumed  to  be  zero  mean,  Gaussian,  and  white 

(uncorrelated  in  time)!  with  !i(t)  an  sxs  positive  semi- 
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definite  syirjnetric  matrix  that  is,  in  general,  piecewise 
continuous  in  t. 

v(ti)  « zero  mean,  white  Gaussian,  measurement  noise  sequence 
Independent  of  w(t)  and  xCt^)  for  all  time  (m-vector). 

The  statistics  of  v(tj^)  are  C v(tj^)3  = 0,  and 
ii(^i)v(tj)]  = R(tj^)  tj^  “ tj 

0 otherwise 

= state  transition  matrix  from  time  tj^.^  time  tj^ 

(lixn  matrix) 

P(t^)  e filter  computer  covariance  matrix  of  state  xC'ti)*  also  of 
the  error  in  the  estimate  of  x(tj^),  prior  to  incorporating 
a measurement  at  time  tj^  (nxn  symmetric  matrix) 

£('t^)  *•  filter  computer  covariance  matrix  of  state  x('ti)»  also  of  I 

the  error  in  the  estimate  of  x(tj^),  after  incorporating  a ^ 

measurement  at  time  tj^  (nxn  symmetric  matrix) 

F(t)  •=  system  dynamics  matrix  (nxn)  or  the  matrix  of  partial 

derivatives  of  fLx(t),tD  with  respect  to  x for  the  extended 
KaLman  filter 

G(t)  “ system  input  matrix  (nxs) 

H(tj^)  «=  system  observation  matrix  at  time  tj_  (mxn)  or  the  matrix 
of  partial  derivatives  of  hLx('t)»t]  with  respect  to  x 
for  the  extended  Kalman  filter 

“ positive  definite  measurement  noise  covariance  matrix  (mxra) 

K(tj^)  = Kalman  gain  matrix  (nxm)  defined  at  tine  tj^ 


Linear  Kalman  Filter  Formulation 

The  linear  Kalman  filter  formulation  presented  In  this  section  Is 
for  a continuous  time  system  model  with  discrete  time  measurement 
updates.  Assume  that  the  system  modeling  has  been  completed  and  that 
the  state  vector  x(t)  satisfies  the  vector  stochastic  differential 
equation: 

x(t)  = F(t)  x(t)  + G(t)  w(t)  (2) 

The  state  equation  is  propagated  forward  in  time  from  the  Initial 
condition  x(^o^‘  Since  the  exact  initial  condition  may  not  be  known, 

f 

it  is  modeled  as  being  a Gaussian  random  variable  with  mean  jCq  and  j 

covariance  P : ^ 

. I 

CsCt,,)]  • io  (3) 

C<s(to)  - - Jo  • w 

The  initial  covariance  matrix  ^ may  be  positive  semldefinlte, 
admitting  exact  knowledge  of  the  initial  conditions  of  some  of  the 
states  or  linear  combinations  thereof.  It  can  be  shown  (Ref  10:157- 
163)  that,  luider  the  assumption  that  x('to)  is  either  deterministic 
or  a Gaussian  random  variable,  the  solution  x(t)  to  lineaur  stochastic 
differential  equations  such  as  Equation  (2)  is  a Gauss-Markov  process, 
l,e,  the  conditional  density  of  x at  time  tj^  based  upon  all  realiza- 
tions of  X through  time  tj^.^  i®  Gaussian  and  completely  determined 

hf  the  process  value  at  t^-i.  Because  the  conditional  density  is 
Gaussian,  It  is  completely  specified  by  its  mean  and  covariance 
(Ref  7:92), 
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Measurements  are  available  at  discrete  tine  points  and  are 
assumed  to  be  of  the  form  of  a linear  combination  of  the  states  and 


corrupted  by  a white  Gaussian  sequence  (Ref  352); 

*(ti)  " H(ti)  x(ti^)  + v(tj^)  (5) 

The  state  estimate  propagates  between  measurements  (from  tine 
ti_i  to  tine  tj^)  according  to: 

and  the  covariance  propagates  according  to; 

P(\)  P(tt.i)  $(ti,t^.3^)^ 

+ ^(tj^.T)  G(t)  dT  (7) 


At  neasiurement  time  t^,  the  estimate  is  updated  according  to 
(Ref  10:233); 

x{tj;)  *=  i(t7)  + K(t^)  [X^  - H(t^)  £(t-)]  (8) 

P(tt)  “ P(ti)  - K(ti)  H(ti)  F(tp  (9) 

where 

K(\)  - P(ti)  iKt^)"^  [H(t^)  P(tp  H(ti)^  + B(t^)T^  (10) 

where  [ 3“^  indicates  the  inverse  of  the  bracketted  matrix  and  Xi 
the  realized  value  of  the  measurement  z(tj^)  at  tine  t^^. 

Under  the  assumption  that  the  adequate  system  model  is  linear, 
and  that  the  dynamic  driving . and  measurement  noises  are  Gaussian  and 
white,  the  Kalman  filter  provides  the  optimal  estimate  x(ti)  of  the 
state  of  the  system  (Ref  10:66,214),  relative  to  many  optimality 
* criteria  with  these  assumptions  x(t^)  is  the  mean,  node,  and  median 
of  the  conditional  density  of  x(tj^),  conditioned  on  the  entire  measure- 
ment history  through  time  t, . The  covariance  of  the  error  conunitted 
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ly  using  x(t^)  as  the  estimate  of  the  state  at  time  tj^  is  denotea 
P(t^),  It  should  be  noted  that  for  a linear  estimation  problem,  the 
covariance  propagation  liquations  (7,9)1,  while  depending  on  the 
sequence  of  H(tj^)  and  R(tj^),  is  independent  of  the  time  history  of 
measurements  ]C2»  •••)•  This  will  no  longer  be  the  case  in  the 

extended  Kalman  filter  formulation. 

The  assumption  that  the  system  can  be  modeled  as  being  driven  by 
white  Gaussian  noise  is  often  well  founded  on  two  accounts.  First, 
it  has  been  found  from  practical  experience  that  the  Gaussian  distri- 
bution provides  a reasonable  approximation  to  observed  random  behavior 
in  certain  physical  systems  (Ref  7*92).  Secondly,  the  central  limit 
theorem  (Ref  7*96)  states  that  if  the  random  phenomenon  that  we  observe 
at  the  macroscopic  level,  is  due  to  the  superposition  of  a large  number 
of  independent  random  variables  on  the  microscopic  level,  then  the 
macroscopic  phenomenon  can  bo  adequately  modeled  as  a Gaussian  random 
variable  (Ref  10s40). 

Extended' Kalman  Filter  Formulation  (Ref  3*179-189) 

The  extended  Kalman  filter  formulation  is  commonly  used  in  esti- 
mation problems  in  which  the  adequate  state  and/or  measurement  equations 
are  nonlinear  rather  than  linear.  Consider,  as  before,  a system  that 
is  continuous  in  time  with  measurements  at  discrete  sampling  times. 
Assume  that  the  system  state  satisfies  the  following  nonlinear  vector 
stochastic  differential  equation* 

*(■*-)  “ TQ£('t),  t]  + G(t)  w(t)  (ll) 
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where  is  a nonlinear  function  of  the  states  and  time  [[in 

general  f could  also  fee  a function  of  deterministic  control  inputs 

where  the  vector  w(t)  of  zero  mean  white  Gaussian  driving 
noises  and  co'/ariance  g(t)  enters  in  a linear  additive  manner.  The 
initial  condition  of  x('to)  modeled  as  a Gaussian  random  variable 
with  mean  ^ and  covariance  P^.  Noise  corrupted  vector  measurements 
of  a nonlinear  function  of  the  states  and  time  are  available  at 
dlsci'ete  times  tj^  as; 

z(tj^)  tj^]  + v(tj^)  (12) 

where  v(t^)  is  a zero  mean  white  Gaussian  sequence  with  covariance 
kernel  (Ref  3*120). 


Cx(ti)  v(tj)] 


H(ti) 

0 


for  1 = j 
otherwise 


(13) 


It  is  assumed  that  the  processes  w and  v are  independent  of  each  other. 
The  filter  propagation  equations  axe; 

i(t/tj^)  =£x(t/t^),  t]  (14) 

i(tiAi)  = i(ti)  (15) 

where  the  notation  x(t/tj^)  means  the  optimal  estimate  of  the  state, 

X,  at  time,  t,  given  the  updated  estimates  up  to  and  including  time 
In  addition,  (the  covariance  is  propagated  approximately  by); 
P(t/tj^)  = £t;i(t/tj^)]  P(t/tj^)  P(t/tj^)  £Et;x(t/ti)]^ 


+ G(t)  a(t)  ^(t)"^  (16) 

P(ti/ti)  = Il(t^)  (17) 

where  F is  the  matrix  of  partial  derivatives  of  f with  respect  to  x, 
evaluated  along  the  tra jectoiy  [[  which  is  propagated  by  means  of  Equation 

(14)]  I 


1 
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The  measurenent  update  is  given  byt 

K(ti)  “ £(^1)  P(tp  ££ti;x(tp] 

+ (19) 

x(t^)  = U\)  + K(ti)CXi  - ]£x(ti),ti]}  (20) 

P(ti)  = P(ti)  - K(t3^)  ££t;i(t;;)]p(t7)  (21) 

vhere 

= i('tiAj^_i)  (22) 

P(ti)  •=  P(ti/tj^_3^)  (23) 

r n ... 

l£tii(tp]  - (24) 

- X “ x(ti) 

The  notation  t^  implies  the  value  of  the  quantity  at  the  instant  prior 
to  the  update  at  time  t^^,  and  t^  is  the  value  of  the  quantity  just  after 
the  update.  The  notation  used  thus  far,  except  as  noted  for  the  ? 
matrix,  is  developed  and  used  in  (Ref  3). 

Unlike  the  conventional  Kalman  filter,  the  extended  Kalman  filter 
gain  and  estimation  error  covariance  matrices  depend  on  the  time  history 
of  x(t/tj^).  Since  a covariance  analysis  could  only  use  an  irominal^^^ 
as  an  approximation  to  x(t/t^)  as  actually  used  by  the  filter  in  its 
online  operation  (the  linearized  Kalman  filter),  the  extended  Kalman 
filter  performance  should  be  verified  by  a Monte  Carlo  simulation. 

In  general,  this  form  of  analysis  is  both  more  time  consuming  and  more 


costly  in  terms  of  computer  usage  than  is  a covariance  analysis. 
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Section  III 


Hontc  Carlo  Simulation  and  Analysis  of  a Sub-Optimal  Kalman  Filter 

Introduction 

As  mentioned  in  Section  I,  a llonte  Carlo  analysis  is  a computer 
simulation  technique  that  can  be  used  to  develop  the  required  (extended) 
Kalman  filter  tuning,  portray  filter  performance  capabilities,  and 
study  the  effects  of  parameter  variations  in  the  true  system  model. 

The  Monte  Carlo  simulation  uses  random  number  generators  and  shaping 
filters  which  generate  random  errors  to  noise  corrupt  both  system 
state  dynamics  and  measurements.  The  noisy  system  measurements  are 
then  processed  by  the  suboptimal  (extended)  Kalman  filter  to  generate 
the  filter  estimates.  If  many  of  these  simulations  are  conducted, 
then  one  can  compute  the  statistics  of  the  error  between  the  filter 
estimate  and  the  truth  model  for  quantities  of  interest. 

Performance  Analysis 

Figure  1 depicts  schematically  a means  of  conducting  the  perform- 
ance analysis  of  a given  (extended)  Kalman  filter  design.  The  truth 
model  by  def'nition  is  the  best,  most  complete  mathematical  model  that 
can  be  developed  to  describe  the  system  under  study.  Such  a truth  model 
is  the  product  of  extensive  study  and  data  analysis  of  the  system. 

As  noted  in  the  figure,  the  truth  model  is  an  n-fState  model,  linear 
or  nonlinear,  driven  by  noise  wt(t)  (assumed  white  and  Gaussian),  that 
generates  the  true  state  values  x^(t).  It  is  assumed  that  the  true 
values  for  the  critical  quantities  of  interest  are  related  to  the  true 
state  values  by  the  vector  function  c^  (a  p vector  - generally  less 


Often,  £^Cx^(t),t]  is  a lin''ar  function  of  x^,  as 

= £(t)x(t)  (26) 


In  fact,  in  many  instances,  the  quantities  of  interest  are  simply  a 
subset  of  the  truth  model  states;  letting  the  components  of  x(t)  be 
ordered  so  that  the  first  p states  are  the  quantities  of  interest, 
this  yields  the  particular  form» 

= Li!^  xt(t)  (27) 

Thus  we  have  access  to  the  "real  world"  values  for  comparison  with  the 
filter  estimates. 


At  discrete  times  (tj^)  the  m-dimensional  measurements* 

+ Xt(ti)  (28) 

are  presented  to  the  Kalman  filter.  The  Kalman  filter  algorithm 
processes  the  measurements  and  produces  the  state  estimates  x(t). 

These  state  estimates  are  related  to  the  estimates  of  critical  quanti- 
ties of  interest  by  the  vector  function  £ (a  p vector): 

i(t)  -c[5(t),t]  (29) 

Thus  the  true  error  committed  by  the  Kalman  filter  (£^)  at  time  t^^, 
before  and  after  measurement  incorporation,  is: 

et(^D  = i(^i)  " (30) 

£t('tt)  = i(tp  - it(^i)  (31) 

If  a feedback  control  system  with  impulsive  correction  is  assumed, 
then  a third  error  is  of  interest  as  wells 

et(tl°)  -it(0  (32) 

where  the  superscript  c denotes  after  the  control  is  applied. 
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The  objective  of  the  performance  analysis  is  to  characterize  the 
error  process  statistically.  This  is  accomplished  in  the  ilonte  Carlo 
simulation  by  generating  many  samples  of  the  error  process  and  then 
computing  the  sample  statistics  directly.  If  enough  samples  are  gene- 
rated, then  the  sample  statistics  should  approximate  the  process 
statistics  very  well  (Ref  , The  sample  statistics  computed 

at  each  point  in  time  are  the  calculated  mean  error  (i),  calculated 
standard  deviation  of  the  error  (s^) , and  the  ensemble  average  of  the 
filter  covariance  diaigonal  terms  • The  calculations  are  made 

over  the  ensemble  of  runs  (n)  for  each  time  point  (tj): 


1 N 


e(tp  = N (33) 

s (t)  = [{  2 [Xi(t  ) - x^i(t  )]^  - U e^(t  (34) 

^ 1=1  J J J 


- i 2^  (35) 

These  equations  can  be  found  in  any  good  statistics  text  (see  for 


exeimple  Ref  11:26,40).  It  is  important  to  note  the  use  of  l/(K-l) 
in  the  expression  for  the  standard  deviation,  Sg,  instead  of  lAs 
this  re ’ ^Its  in  an  unbiased  estimator  of  variance  before  the  [ 
is  taken  and  is  customarily  used  for  sampled  data  (for  small  or 
moderate  N), 

One  useful  output  of  the  Monte  Carlo  analysis  is  a plot  of  the 


filter  estimate  of  the  standard  deviation,  t . ) . along  with  the 

corresponding  computed  standard  deviation,  Sg(*j)  '*'j 

• interest.  The  Kalman  filter  is  tuned  by  minimizing  the  computed  error 
standard  deviations.  As  a rule  of  thumb,  good  tuning  is  achieved 
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when  the  filter  estimated  standard  deviations  are  approximately  equal 
to  the  computed  standard  deviations.  If  the  problem  imder  investiga- 
tion is  sensitive  to  small  changes  in  the  filter  tuning  parameters 
(]p(t^),  and  time  histories  of  ^(t)  and  B(t)3  then  the  investigator 
must: 

1,  Conduct  the  Monte  Carlo  simulation 

2,  Examine  the  standard  deviation  plots 

3,  Determine  tuning  parameter(s)  to  be  changed 

4,  Repeat  1 through  3 until  adequate  tuning  is  achieved 

Thus  it  can  be  readily  visualized  that  precise  tuning  of  the  Kalman 
filter  by  Monte  Carlo  techniques  is  very  costly  in  terms  of  effort  and 
time. 

A second  useful  output  of  the  Monte  Carlo  analysis  is  a plot  of 
the  mean  error  versus  time.  This  plot  is  usually  obtained  with  one 
standard  deviation  bounds  plotted  above  and  below  the  mean  error. 

This  plot  is  useful  to  determine  if  the  suboptimal  Kalman  filter 
prov'ides  the  desired  accuracy.  If  the  filter  starts  from  initial 
conditions  other  than  the  truth  model  initial  conditions,  then  the 
filter  performance  to  non-zero  initial  error  can  be  observed.  Bias 
errors,  a problem  of  significance  to  many  extended  Kalman  filter 
designs,  also  become  readily  apparent  from  the  plots  of  this  form. 

Digital  Simulation  of  the  Truth  Model  and  Filter  Model 

In  order  to  simulate  the  system  on  a digital  computer,  the  problem 
• must  be  divided  into  segments  which,  when  processed  sequentially  over 
a given  time  period,  will  provide  a realistic  portrayal  of  filter  use 


i 

i 


j 

i 
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in  actual  on-line  inpleracntation,  A convenient  time  period  is  the 
sample  period  between  tines  at  which  neasurenents  become  available.  J 

I 

Then  the  entire  procedure  is  iterated  over  some  desired  time  interval  5 

of  interest.  The  segments  are;  j 

■! 

1.  Time  propagation  of  system  state  equations  | 


2.  Tine  propagation  of.  filter  state  equations  | 

3.  Time  propagation  of  filter  covariance  equations  ^ 

4.  Update  of  filter  covariance  matrix  ^ | 

i 

5.  Generation  of  measurements  for  the  filter  ] 

] 

6.  Update  of  filter  states  i 

7.  Application  of  impulsive  feedback  to  "update"  system  states  ] 

1 

8.  Reset  the  filter  states  after  the  feedback  | 

The  first  segment  is  simulated  by  integrating  the  n^-dimensional  ! 

system  (truth  model)  stochastic  differential  equation;  1 

= f^[x.,.(t),t]  + G^(t)w^(t)  (36)  l| 

j 

from  time  tj^  to  time  of  the  next  measurement  tine  The  Integra-  I 

tion  is  typically  performed  by  either  Euler  or  Runga-Kutta  numerical  j 

j 

integration  methods  with  a specified  Integration  step  size.  Integrating  ; 

the  deterministic  part  of  the  differential  equation  is  readily  apparent,  • j 

however,  the  additive  white  Gaussian  noise  w.(.  is  not.  The  following  j 

derivation  shows  how  the  solution  to  the  stochastic  differential  | 

equation  can  be  approximated.  ; 

To  provide  insight  into  the  approximate  solution  form,  consider  | 

1 

first  a linear  stochastic  differential  equation  of  the  form;  | 


I 


vhere  w(t)  is  a zero  mean  white  Gaussian  noise  with  arw(t)wCt+t)'^~]  ' 

^(t)6(t).  An  equivalent  discrete  time  system  equation  would  bej 

^i+1 

x(ti+i)  = i(t,+i,t.)  x(t.)  + / t(ti+i,t.)  G(t)h(t)  dT 

(38: 

= Vi'\^  ^^^i^  (39: 

where  satisfies  both  of  the  following  relationships: 

^^(t,ti)  = F(t)  S(t,tj^)  for  t in  the  interval 

= I (41) 

For  this  discrete  time  system,  a zero  mean  white  Gaussian 

discrete- time  noise  with  strength: 

^ G(t)  £(t)  £(1')"^ 


^(■ti+3^,T)  dT 


2d^"i) 


] = 0 for  tj^  tj  (44) 

Now  if  the  integration  step  size  is  small  relative  to  system  character- 
istic response  tine  and  constant  (At  sec  in  duration)  then  the  following 
approxinstions  can  be  made; 

i(ti  i.t^)  = I + F(t^)  At  (45) 

i G a dT  = CG(t  ) £(ti)  G(t At  (46) 

2 

where  all  terms  of  order  At  or  higher  have  been  neglected.  Thus  one 
approximate  simulation  of  the  continuous  tine  system  would  be: 

x(Vi)  = [I  + F(ti)  At]  x(t^)  + w^(t^)  (47) 

•=  + £(^1)  x(\)  dt  + 2iji(tj^)  (48) 


where  is  described  by: 

“ LG(ii)  9S^\)  5.(ii)^]  (49) 

Thus  the  differential  equation  is  integrated  using  Equation  (46) 
and  noise  generators  utilized  for  the  i[<j(tj^)  tern.  Higher  oitier  inte- 
gration tecbjiiques  can  be  applied  to  the  nonhomogeneous  portion  of  the 
differential  equation,  with  the  same  >i{j^(tj^)  as  above  used  as  the  forcing 
terms.  It  should  be  noted  that  this  technique  is  directly  extendable 
to  nonlinear  equations  of  the  form: 

x(t'^  = f[x(t),t]  + G(t)  w(t)  (50) 

Segment  two  is  performed  by  calculating  the  filter  equation  (6). 

The  filter  covariance  propagation  Equation  (?)  is  calculated  to 
accomplish  segment  three.  Segment  four  is  performed  by  solving 
Equations  (lO  and  9).  Segment  five  uses  Eqi;iation  (5)  to  generate  the 
m-dimensional  measurements.  Noise  generators  are  utilized  to  provide 
the  v^(tj^)  term.  The  filter  update  Equation  (8)  is  calculated  to 
perform  segment  six. 

A Kalman  filter  will  often  be  implemented  in  an  indirect  feedback 
configuration,  in  which  filter  error  state  estimates  are  fed  back  to 
the  actual  system  to  correct  the  actual  er.  :r.  If  the  system  correc- 
tion can  be  performed  rapidly  enough  compared  to  the  filter  update 
period,  then  the  feedback  can  be  simulated  as  an  impulsive  change. 

After  application  of  the  correction,  assuming  linear  feedback,  the 
truth  model  state  process  becomes  (Ref  3*6-76): 
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The  filter  should  be  "told"  that  this  feedback  to  the  system  has 


occurred,  so  its  state  estimate  is  modified  as; 

=^{t^)  - D(tj^)  x(t^)  (52) 

= [I  - £(ti)]  x(ti)  (53) 

Note  that  if  the  filter  is  estimating  error  states  then  the  D(tj^) 
matrix  will  typically  be  an  identity  matrix.  Segments  seven  and 
eight  axe  performed  by  calculating  the  results  of  Equations  (51) 
and  (52). 

Each  of  these  segments  are  performed  sequentially  at  the  measure- 
ment update  time  and  the  simulation  time,  system  state,  filter  state, 
and  filter  covariance  values  saved  on  a data  tape.  This  process  con- 
tinues until  the  simulation  final  time  is  reached,  thus  completing 
one  run  of  a Monte  Caxlo  simulation.  After  many  runs  of  the  Honte 
Carlo  simulation  have  been  pei-formed,  then  the  data  tape  can  be  read, 
the  sample  statistics  calculated  using  Equations  (33“35)  and  computer 
plots  created. 


SECTION  IV 


DESCRIPTION  OF  THE  GENERALIZED  MONTE  CARLO  ANALYSIS  PROGRAll  (MCAP) 
Introduction 

The  general  Monte  Carlo  Analysis  Program  (flCAP)  is  a modular 
program  designed  to  perform  a Monte  Carlo  performance  analysis  of  a 
given  (extended)  Kalman  filter  design.  The  program  is  composed  of 
the  main  program,  problem  independent  (constant)  subroutines,  and 
a minimum  of  8 user  defined  problem  dependent  (variable)  subroutines. 

The  main  program  acts  as  an  executive  to  sequence  the  simulation 
through  the  performance  analysis.  The  MCAP  subroutines  can  be  divided 
into  three  major  categories:  l)  data  input  subroutines,  2)  data  output 
subroutines,  and  3)  "workhorse"  subroutines.  The  data  input  subroutines 
read  in  the  problera  dependent  parameters  and  control  parameters  for 
MCAP  execution.  The  output  from  MCAP  consists  of  data  printed  on  the 
computer  line  printer  and  a data  file  (Tape  3)  containing  the  sample 
realizations  for  the  many  simulation  runs.  The  amoimt  of  printed 
output  is  determined  by  the  control  parauneters  specified  in  the  data 
input.  The  output  data  file  (Tape  3)  is  used  by  the  statistics/plot 
(STATPLOT)  program  to  calculate  the  statistics  of  the  true  error  noted 
in  Section  III  (Equations  33  through  35)  and  to  prepare  a plot  file 
for  the  Calcomp  (or  other)  plotting  device.  The  "workhorse"  subroutines 
are  executed  in  the  sequence  managed  by  the  main  program  to  perform 
many  functions : 

1,  Initialize  the  Monte  Carlo  performance  analysis 

2,  Initialize  the  simulation  run 

3,  Propagate  the  system  state  equations 
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4.  Propagate  the  filter  state  equations 

5.  Propagate  the  filter  covariance  equations 

6,  Update  the  filter  covariance  matrix 

7,  Generate  system  meas\u:ements  for  the  filter 

8,  Update  the  filter  states 

9.  Apply  feedback  (if  applicable  to  the  problem)  to  reset  the 

system  and  filter  states 

KCAP  vas  developed  using  the  General  Covariance  Analysis  Program 
(GCAP)  by  the  Air  Force  Avionics  Laboratory  (Ref  2)  as  a baseline. 
Consequently,  the  goals  of  having  MCAP  similar  in  structure  to  GCAP 
and  having  (some)  problem  dependent  subroutines  applicable  to  both 
GCAP  and  fICAP  are  met.  Because  MCAP  is  composed,  to  a large  extent, 
of  modified  GCAP  subroutines;  users  already  familiar  with  GCAP  should 
not  have  any  difficulty  in  applying  hCAP  to  their  problem. 

As  noted  in  Section  I,  in  general,  the  covariance  analysis  is 
limited  because  stringent  assumptions  are  necessary  for  the  analysis 
results  to  be  a valid  depiction  of  the  error  characteristics.  In 
particular  GCAP  assumes  that: 

1.  The  filter  auid  system  are  modeled  by  a set  of  linear  differen- 

tial equations. 

2.  The  measurements  used  by  the  filter  are  a linear  combination 

of  the  system  states  with  additive  noise-corrupted  measure- 
ments , 

3.  The  states  being  modeled  in  the  filter  are  error  states. 


J 
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4,  Total  impulsive  feedback  control  is  applied  to  the  system 
(error)  states. 

These  stringent  assumptions  limit  the  application  of  GCAP  to  carefully 
posed  problems.  None  of  these  assumptions  are  necessary  for  problems 
to  be  investigated  using  MCAP,  Thus,  at  the  expense  of  additional 
computer  time,  MCAP  can  be  used  to  solve  a larger  class  of  estimation 
problems, 

MCAP  Program  Structure 

The  performance  analysis  of  the  Kalman  filter  is  conducted  under 
the  manaigement  of  the  main  executive  program.  Figure  2 illustrates 
the  structure  of  flCAP.  There  are  generally  three  levels  of  subrout j.nes 
within  this  structure.  The  first  level  subroutines  are  coded  on  the 
figure  to  note  the  main  functions  performed  and  if  the  subroutine  is 
problem  dependent  (a  variable  subroutine  to  be  specified  by  the  user) . 
The  function  code  is  located  below  the  subroutine  name.  If  the  sub- 
routine is  problem  dependent  then  the  subroutine  name  is  underlined 
(for  example  XSDQT).  The  coding  used  on  this  figure  isi 
I - Data  input  function 
0 - Data  output  function 

W1  - Initialize  the  Monte  Carlo  performance  analysis 

W2  - Initialize  the  simulation  run 

W3  - Propagate  the  system  state  equations 

W4  - Propagate  the  filter  state  equations 

W5  ••  Propagate  the  filter  covariance  equations 

W6  - Update  the  filter  covariance  matrix 


W7  - Generate  system  measurements 

2>’+ 


Figure  2.  ]ICAP  Program  Structure  (con't) 
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W8  - Update  the  filter  states 


W9  - Apply  feedback  to  update  the  filter  and  system  states 
Several  of  the  subroutine  names  are  repeated  to  Illustrate  the  program 
structure  and  sequence  of  subroutine  execution  better. 

The  first  task  WOAP  perforins  is  to  read  in  the  problem  control 
parameters,  problem  dependent  parameters,  and  assign  memory  storage 
locations.  Function  W1  is  than  psrfcmed  in  order  to  set  the  print 
and  plot  control  "switches"  and  initialize  the  random  number  generator. 
The  simulation  run  is  initialized  (^2)  by  calculating  the  integration 

constants  and  initializing  the  program  time  parameters. 

Ij 

I One  run  of  the  llonte  Carlo  simulation  is  performed  by  sequencing 

through  a cycle  of  functions  W3  through  W9  at  a rate  equal  to  the 
measurement  interval  until  the  final  problem  time  is  reached.  Through- 
out this  run,  the  data  output  fvinction  (o)  is  perfonp.ed  to  print  the 
data  selected  by  the  user  with  the  problejTi  control  parameters. 
Subsequent  runs  are  performed  by  repeating  functions  W2  through  W9 
until  the  user  specified  number  of  runs  is  completed. 

The  rest  of  this  section  expands  on  this  brief  description  and 
provides  a description  of  the  variable  subroutines  required  for  iiCAP, 

A good  description  of  the  constant  subroutines  is  contiiined  in  the 
GCAP  users  guide  (Ref  2)  and  is  not  repeated  here. 

j Data  Input  Function  (l).  The  variables  listed  and  described  in 

i|  Table  I are  the  input  parameters  required  for  MGAP.  The  normal  input 

!! 

j media  for  the  problem  variables  is  the  caird  deck  structure  shown  in 

;[  Figure  3«  However,  the  user  may  elect  to  read  the  values  of  I’^', 
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Table  I 


Definitions  and  Format  of  iiCAP  Input  Parameters 


Parameter 

Description 

Format 

(title) 

Problem  title  (limited  to  80  characters) 

6a10 

$USRCTL 

Problem  dependent  control  parameters 
are  contained  in  this  namelist 

Namelist 

Nl, 

Suboptimal  filter  dimension 

Namelist 

N2, 

Truth  model  ('system  dimension) 

Namelist 

M, 

Number  of  measurements 

Namelist 

Initial  problem  time  (seconds) 

Namelist 

TF, 

Final  problem  time  (seconds) 

Namelist 

DTUP, 

Measurement  update  interval  (seconds) 

Namelist 

DTINT, 

Integration  Interval  (seconds) 

Namelist 

IPASS, 

Number  of  runs  in  Monte  Carlo 
simulation 

Namelist 

I3E£D, 

Initial  seed  for  random  number 
generation 

Namelist 

IRNPRCT, 

Number  of  runs  per  printed  output 

Namelist 

IPRCT, 

Number  of  integration  intervals  per 
printed  output 

Namelist 

IPRSYH, 

if  0 - do  not  print  lower  triangular 
portion  of  syruiietric  matrix  PF 
if  1 - print  lower  diagonal  of  symmet- 
ric matrix  FF 

Namelist 

IPRK, 

if  0 - do  not  print  Kalman  gain  matrix 

K 

if  1 - print  Kalman  gain  matrix  K 

Namelist 

IPLCT, 

Nvunber  of  integrations  intervals  per 
plotted  output  (must  be  1 to  get 
plot) 

Namelist 

Table  I (cont'd) 


IPFILE, 

1 

if  0 - usual  data  input  from  cards 
if  1 - input  restart  values 

TIMUD{/,  XF,  XS,  and  PF  from  Tape  11 

Namelist 

TIMUD5!f, 

Starting  value  of  the  update  time  in 
seconds  — imused  if  IPFILE  = 1 

Namelist 

TIJiNTgf, 

Starting  value  of  the  integration 
time  in  seconds  — unused  if 

IPFILE  = 1 

Namelist 

$ 

End  of  Namelist 

Namelist 

NZFF 

1 

Number  of  non-sero  terms  in  the 
filter  F matrix 

13 

INDFF 

The  coordinate  pairs  of  the  non- 
zero terms  in  the  filter  F matrix 

8(2I3,4X) 

KZQF 

Number  of  non-zero  terms  in  the 
filter  ^ matrix 

13 

INCqF 

The  coordinate  pairs  of  the  non-zero 
terms  in  the  filter  ^ matrix 

8(213,^0 

NZIiF 

Number  of  non-zero  terms  in  the  filter 

H matrix 

13 

INDHF 

The  coordinate  pairs  of  the  non-zero 
terms  in  the  filter  H matrix 

8(2I3,4X) 

NZRF 

Number  of  non-zero  terms  in  the  filter 

R matrix 

13 

INDRF 

The  coordinate  pairs  of  the  non-zero 
texviiii  in  the  fix vez  ^ watr.  ‘c 

8(213,4/:) 

XFijS 

if  IPFILE  = 0,  each  filter  state  ini- 
tial condition  value  (N1  in  niMber) 

E20.0 

if  IPFILE  = 0,  each  system  state  ini- 
tial condition  value  (N2  in  number) 

E20.0 

FF^C 

if  IPi'TLE  = 0,  each  filter  covariance 
diagonal  element  initial  condition 
value  (N1  in  number) 

i 

ElO.O 
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TIMUJJ^,  XF{2f,  and  from  the  file  neimed  "Tape  11"  as  shown  in 

Figure  4.  This  option  facilitates  "restarting"  the  simulation  at  the 
simulation  time  where  a previous  execution  was  terminated. 

Concurrent  with  the  data  input,  the  variables  N2  and  li  undergo  a 
"reasonableness"  test  in  order  to  insure  program  memory  is  not  exceeded. 
If  the  values  of  these  variables  pass  the  test  criteria  (M  less  than 
21  and  HZ  less  than  40l) , then  J-'CAP  begins  to  assign  storage  locations 
within  the  unlabeled  CCIIIION  block  A.  These  storage  locations  are  used 
to  store  the  indices  of  the  filter  state  matrices  (F,  K,  and 
When  this  process  is  completed  for  the  filter  state  matrices,  "CAP 
generates  the  memory  location  pointers  for  the  first  word  locations 
of  the  initial  values  of  the  filter  states,  the  system  states,  and 
the  covariance  matrix;  the  time  varying  values  of  the  filter  states, 
the  system  states,  and  the  covariance  matrix;  the  time  varying  values 
of  the  filter  state  derivatives,  the  system  state  derivatives,  and 
the  filter  covariance  derivatives;  the  filter  state  matrices  (F, 

H,  and  R),  the  Kalman  gain  matrix;  and  memory  locations  for  temporary 
intermediate  calculations.  This  elaborate  memory  assignment  scheme 
enables  the  reuse  of  storage  locations  while  the  simulation  is  executing. 
However,  th.>  user  should  be  cautioned,  that  this  "feature"  may  cause 
problems  if  values  are  needed  in  the  user-supplied  vsiriable  subroutines. 
The  user  then  has  the  option  of  over-riding  the  HCAP  storage  assignment 
[by  modifying  the  memory  assignment  subroutine  (STRVAL)]  or  storing 
the  values  of  the  variables  in  the  variable  subroutines  as  labeled 
COraON. 
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1 

Initialize  the  Performance  Analysis  . Once  the  storage 
allocations  are  made,  MGAP  sets  the  print  and  plot  control  "switches" 
according  to  the  problem  control  parameters  (IRIJFRGT,  IPRCT,  IFRoYil, 

IPRX,  and  IPLGT).  The  random  number  generator  "seed"  is  initialized 

to  the  value  specified  by  the  user  (ISSED).  This  seeding  of  the  random 

number  generator  is  for  Ilonte  Garlo  simulation  repeatability  from  one 

set  of  runs  to  another.  Next,  the  user  supplied  subroutine  UihIN  is  ^ 

called.  The  subroutine  USRIN  may  be  employed  in  a variety  of  ways, 

USRIH  may  be  used  to  input  parameters  which  change  between  honte  Garlo  i | 

simulations  but  which  remain  unchanged  within  a single  simulation.  Or 
the  user  may  employ  USRIN  to  create  special  initial  conditions,  such  as 
a nondiagonal  FF^^,  Whatever  the  usaige,  USRIN  is  a convenient  method 
which  permits  a user  to  customize  MGAP  input  data  to  the  specific 
pi-oblem. 

Initialize  the  Simulation  Run  (W2) , MGAP  then  initializes  the 
time  dependent  values  for  the  filter  states,  the  system  states,  and 
the  covariance  matrix  to  the  values  specified  in  the  data  input.  The 
two  user  supplied  subroutines  XFDOT  and  X3D0T  are  called  to  allow  these 
initial  values  to  be  modified  with  each  simulation  run.  The  integra- 
tion constants  and  the  program  time  parameters  are  then  initialized, 

MGAP  is  now  reauiy  for  a simulation  run. 

Propagate  the  Derivative-  Equations  (W3  through  W5),  One  run  of 
the  Monte  Garlo  simulation  involves  sequencing  through  a cycle  of 
functions  W3  through  W9  at  a rate  equal  to  the  measurement  update 

242 


A 


Interval}  until  the  final  problem  time  is  reached.  Functions  W3 
throu/jh  W5  require  the  extrapolation  of  the  system  state  values,  the 
filter  state  values,  and  the  covariance  matrix  from  the  current  simu- 
lation time,  TIMS,  (which  hcas  an  initial  value  of  1^)  until  the  first 
update  time,  TIliUD.  TIliUD  is  the  sam  of  TIMS  and  the  measurement  up- 
date interval,  DTUP.  Subroutine  IliTijCiX  pi’opagates,  in  turn,  the  system 
state  values  and  the  filter  state  values  between  TIMS  and  TIMUD, 

INTBSX  accomplishes  this  propagation  by  successively  integrating  the 
system  and  filter  propagation  equations  [Equations  (6)  or  (l4)],  where 
each  integration  is  DTINT  long.  The  integration  is  performed  using  a 
fourth  order  Runge-K  tta  numerical  integration  method.  In  the  case 
of  the  system  state  equations,  the  zero  mean  white  Gaussian  noise 
sequence,  Wjj^(tj^)  [Equation  (48)]]i  is  added  to  the  system  state  values 
at  this  DTINT  rate.  Subroutine  INTEGP  propagates,  in  a similar  manner, 
the  covariance  matrix  from  TIMS  to  TIMUD, 

In  order  to  perform  the  fourth  order  Runge-Kutta  calculations,  both 
INTBGX  and  INTZCP  require  values  for  the  derivative  equations  at  the 
following  time  points:  TIME,  TIME+DTINT/2,  TIME+DTIRT/2,  and 
TIflE+DTINT.  (Note  that  the  fo\arth  order  Runge-Kutta  method  requires 
a derivative  value  twice  at  the  midpoint  of  the  integration  interval). 
For  the  system  and  filter  derivative  states  this  is  accomplished  by 
calling  either  subroutine  X3D0T  or  XFDOT.  The  calculations  for  f , on 
the  other  hand,  requires  values  for  the  F and  ^ matrices,  vihich  in  turn, 
derive  their  values  from  (possibly)  time  vaaying  parameters.  These 
time  varj'ing  parameters  are  obtained  when  INTi)XJP  calls  the  user  sub- 
routine Tl'AJ  at  each  tine  point  noted  above.  The  time  varying 
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pairameters  are  transmitted  between  TRAJ  and  user  subroutine  FL'rilAT 
via  a user  defined  labeled  COIilION.  After  the  call  to  TRAJ,  IUTjJP 
calls  f 'LTl lAT  two  separate  times  in  order  to  obtain  the  current  (non- 
zero) values  of  the  F and  ^ matrices  respectively.  IllTiiJP  then  cal- 
culates P and  extrapolates  the  covariance  matrix  to  the  next  time 
point,  using  Equation  (?)  or  (l6). 

Update  the  Filter  Covariance  Matrix  (V?6) . When  the  TIMUD  time  is 
reached,  function  can  be  performed.  In  order  to  update  the  filter 
covariance  matrix,  the  Kalman  gain  matrix,  K,  must  be  computed  (see 
Section  II).  However,  this  requires  the  current  values  of  the  filter 
H and  R matrices.  These  values  are  obtained  when  subroutine  UPUATFP 
calls  FLTMAT  tvo  separate  times.  Subroutine  UPDATFP  then  calls  other 
subroutines  in  order  to  calculate  K [Eq’jation  (lO)  or  (19)3  ^P" 

date  the  filter  covariance  matrix,  P [^Equation  (9)  or  (21)]. 

Generate  the  System  Measurements  (VJ?) . The  system  measurements  anre 
generated  by  the  user  supplied  subroutine  MZAS.  MSAS  calculates  the 
measurements  using  the  current  system  state  values  corrupted  by 
additive  noise  [^Equation  (5)  or  (12)]. 

Update  the  Filter  and  System  States  (VF8  and  W9).  The  filter 
estimates  (XF^-)  after  the  measurements  obtained  are  generated  by  the 
user  supplied  subroutine  UPDATXF . UPDATXF  generates  XFt  by  performing 
the  calculations  indicated  in  Equation  (20)  [see  Section  II ],  This 
routine  is  user  supplied  in  order  to  accomodate  the  non-linear  vector 
fraction  h.  After  the  XFh-  estimates  have  been  obtained,  the  user 
supplied  subroutine  CQBH5CT  is  called  to  apply  Impulsive  feedback  to 


(selected)  system  states.  (Note  that  if  no  feedback  correction  is 
required  for  the  problem  under  study,  subroutine  CORRxjGT  will  simply 
have  a RETURl'I  and  £2JD  statement) . The  feedback  is  applied  by  modi- 
fying the  system  state  values  using  Equation  (51)  and  the  filter  state 

values  with  Equation  (52).  ! 

! 

Additional  Simulation  Runs.  Throughout  the  simulation  run,  data  I 

is  printed  (values  of  PF,  K,  and  Z)  according  to  the  user 

specified  values  of  the  printer  control  parameters  (IHNPRCT,  IFRCT, 

IPRSYM,  and  IPRK).  In  addition,  at  the  end  of  a simulation  run  (if 
IRNPRGT  is  greater  than  the  current  run  number)  MCAP  prints  the  lower 
triangular  values  of  the  final  ^ matrix  and  "RUN  KUilRER  xx  COMPLETE". 

If  additional  simulation  runs  are  required,  then  functions  W2  through 
W9  are  repeated  until  the  specified  number  of  runs  have  been  completed. 

Variable  Subroutines 

As  mentioned  earlier,  hCAP  requires  a minimum  of  8 problem 
dependent  (variable)  subroutines  to  perform  the  Monte  Carlo  simulation. 

These  variable  subroutines  are: 

1.  XFDOT 

2.  XSDOx 

3.  THAJ 

4.  FLTi’AT 

5.  HEAS 

6.  UPDATXF 

7.  CORRECT 

8.  USRIN 
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These  subroutines  are  designed  to  meet  most  problem  requirements  and 
to  be  used  easily  by  users  having  only  a basic  knowledge  of  Fortran 
programming. 

The  essential  elements  of  each  variable  subroutine  are  described 
below.  In  order  to  illustrate  the  usatge  of  these  subroutines,  the 
next  section  presents  an  example  using  a linear  Kalman  filter  which 
should  clarify  this  discussion. 

XFDOT.  Subroutine  XFSCT.  major  function  is  to  provide  I-iCAP  with 
the  derivative  values  of  the  filter  states  at  each  time  point  of  the 
Runge-Kutta  integration.  The  SUBROUTINE  statement  contains  the 
following  arguments  (in  the  order  shown):  X,  XUOT,  Nl,  INIT,  KOIS, 

DT,  and  TIIIE,  X contains  the  current  values  (an  Kl-dimensioned  vector) 
of  the  filter  states.  The  calculated  values  of  the  filter  state  deri- 
vatives are  returned  to  the  calling  program  (INTEGX)  in  XiX)T  (an  Nl- 
dimensioned  vector),  Nl  is  the  dimension  of  the  filter.  INIT  and 
NOIS  are  declared  logical  variables  to  select  one  of  the  three  sec- 
tions of  this  subroutine  code.  INIT  is  true  only  if  the  subroutine 
has  been  called  during  function  W2,  NOIS  is  false  for  each  time  point 
of  the  Runge-Kutta  integration.  At  the  end  of  the  integration  inter- 
val, The  current  value  of  the  simulation  time  is  contained  in  TIME. 

The  first  section  (when  INIT  is  true)  of  iLFiXDT  is  called  once 
per  simulation  run  (by  INTiOl)  and  is  used  to  initialise  constants, 
unique  to  this  subroutine.  The  second  section  (where  INIT  and  NOIS 
are  false)  is  used  to  calculate  the  values  for  the  derivatives  of  the 
filter  states  using  either  x^  = £(2if»'t)  "the  nonlinear  filter  or 
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^ = F 2^  for  the  linear  filter.  The  third  section  is  noimally  net 
used,  however,  the  option  remains  for  the  user  to  modify  filter 
state  values  at  the  end  of  the  integration  inter\'al. 

Table  II  lists  the  req^uired  Fortran  statements  for  subroutine 
XFDOTa 

Table  II 

Required  Statements  for  XFDOT 


SUBROUTINE  XFDOT  (X,XD0T,N,IKIT,K0IS,DTINT,TIiI3) 
DIJIEIiSION  X(l),XDOT(l) 

LOGICAL  INIT,N0I3 
IFj.NOT.  INIT)  GO  TO  200 
1_ subroutine  constants] 

200  CONTINUE 

IF_(N0I3)  GO  TO  300 

[_derivative  calculations] 

RETURN 
300  CONTINUE 

[user  option] 

RETURN 

END 


X3D0T.  Subroutine  X3D0T  provides  IICAP  with  values  for  the  deriv- 
ative of  the  system  states  (of  N2  dimension)  and,  at  the  end  of  the 
integration  interval,  adds  the  zero  mean  white  Gaussian  noise,  to 
the  system  state  values.  Subroutine  XSDOT  is  identical  in  structure 
to  XFDOT  described  above.  The  only  change  to  Table  II  is  to  modify 
the  subroutine  name  and  replace  the  user  option  section  with  the  calcu- 
lations necessary  to  evaluate  w^  and  add  it  to  the  system  states.  The 
calculations  performed  in  the  second  section  again  uses  either  the 
nonlinear  equation  Xg  = fCxgft)  or  the  linear  equation  x^  = F^  2Ss* 

The  noise  sequence  w^(t^)  covariance  is  found  by  Equation  (^9). 
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The  vector  is  calculated  using  the  IICAF  constant  subroutine 
NOISE,  Subroutine  NOISE  returns  a Gaussian  random  number  with  mean 
and  standard  deviation  specified  on  each  subroutine  call.  The  argu- 
ments for  NOISE  arej  KIS,  XKEAN,  and  XNOISS.  Where  RIIS  is  the 
standard  deviation  of  the  desired  random  number,  XMEAN  is  the  desired 
mean  of  the  random  number,  and  XN'OISE  is  the  returned  random  number. 

If  the  system  dynamics  noise,  w,  elements  are  independent  of  each 
other,  then  matrix  ^ of  Equation  (49)  is  diagonal.  Also,  often 
G is  the  identity  matrix  so  the  term  in  the  brackets  can  be 

written  or  simply  reidentified  The  standard  deviation  for  the 
random  number  is  simply  the  square  root  of  the  ^kk  element  times  the 
square  root  of  At  (DTINT).  Table  III  shows  typical  Fortran  statements 
ij  which  calculate  a-nd  add  it  to  the  system  states  (for  zero  mean 

i independent  white  Gaussian  noise  terms) . If  nondiagonal  G'^G  is 

‘I 

i desired,  this  can  be  handled  bj'"  means  of  transformations  of  indepen- 

|| 

;i  dent  noises  using  square  root  matrices  [^see  Ref  (3s7“6o)3. 


I 

i 


Table  III 

Calculation  of  w^j  for  the  System  States 


2 

for  noise  with  Q..  --a 

XMEAN  = 0. 

RH3  ---  SIGNiA  ^ SQRT(  DTINT) 

CALL  NOISE  (RMS,  XMEAI^,  XNOISE) 

X(J)  = X(J)  + XNOISE 

for  time  correlated  noise  with  = 0 /2t 

where  x(k)  is  a shaping  filter  state; 

XHEAI^  = 0. 

RhS  = SIGMA  * S«iKT(DTINT/(2.*TAU)) 

CALL  NOISE  (RilS,  XMEAN,  Xiv'OISE) 

X(K)  ==  x(k)  + :<I.0ISE 
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TRAJ,  Subroutine  TiiAJ  provides  trajectory  constants  or  time 
varying  parameters  required  by  FLTilAT.  The  requirements  of  TRAJ 

for  MCAP  are  the  same  as  those  for  GCAP,  The  only  argument  for  TIIAJ  i 

is  StlTCON.  SiTCON  is  a declared  logical  variable  to  select  either 

the  tine  invariant  (SETCON  true)  or  the  time  varying  (SETGON  false) 

section  of  the  subroutine.  SETGON  is  true  at  the  beginning  of  each 

integration  interval  and  TRAJ  ca,lcula.tes  the  time  invariant  parameters. 

For  the  rest  of  the  1 itegration  interval  TRAJ  calculates  the  time 
vary'ing  parameters.  Instead  of  computing  the  trajectory  dependent 
parameters  as  the  simulation  is  progressing,  the  user  may  desire  to 
read  the  trajectory  parameters  from  either  a magnetic  tape  or  disk 


unit.  If  this  prestored  parameters  method  is  utilized,  then  the  user  | 

must  Insure  that  TRAJ  transfers  data  to  FLTilAT  at  a rate  which  corre-  j 

I spends  to  twice  the  value  of  DTINT.  Table  IV  lists  the  required  j 

j Fortran  statements  for  subroutine  TRAJ . < 

;j  Table  IV 

I Required  Statements  for  TRAJ 


SUBROUTINE  TRAJ  (SETGOK) 

LOGICAL  Si,TGON 

IF  (.NOT.  SETCON)  GO  TO  100 

[time  invariant  calculations] 
100  CONTINUE 

[time  varying  calculations] 
RETURN 
END 


FLTi'lAT.  oubroutine  FLT?'AT  provides  MGAP  with  nonzero  values  of 


the  filter  watrices  F,  H,  and  R in  turn.  The  F matrix  is  calculated 
using  Equation  (l8)  for  the  nonlinear  filter  or  simply  F in  the  linear 
filter,  is  found  by  calculating  the  term  in  the  brackets  of 
Equation  (46).  The  H matrix  is  calculated  using  Equation  (24)  for 
the  nonlinear  filter  or  simply  H in  the  linear  filter.  R is  the 
measurement  noise  covariance  matrix.  The  requirements  of  FLT;iAT  for 
ilCAP  are  identical  to  those  of  GGAP.  The  subroutine  a,rguments  for 
FLTiAT  are:  SETGON,  lATRO,  A,  and  HID.  3ETC0K  is  again  a declared 
logical  variable  used  to  select  either  time  invariant  or  tine  varying 
calculations  for  the  four  matrix  elements.  The  integer  variable  iATIIO 
determines  the  block  of  coding  accessed  to  calculate  the  nonzero 
values  of  the  desired  matrix.  The  subroutine  returns  the  nonzero 
values  in  the  variably  dimensioned  vector  A.  The  vector  IKD  contains 
the  packed  row-colum.n  pairs  of  indices  which  corresponds  to  the  matrix 
elements  retvjrncd  in  the  A array.  Table  V lists  the  required  Fortran 
statements  for  subroutine  FLTliAT. 


Table  V 

Required  Statements  for  FLTJiAT 


SU3R01JTINS  FLTIAT  (SETGON,  HATNO,  A,  IND) 
DIMEKSION  A(1),  IND(1) 

LOGICAL  3ETC0K 

GO  TO  (101,  201,  301,  401),  KATRO 
101  COliTINUE 

IF  (.IjOT.  33TC0N)  GO  TO  110 

[_time  invariant  calculations  of  F] 
110  CONTINUE 

Ltime  varying  calculations  of  F] 
RETURN 
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Table  V (cont'd) 


201  CONTINUE 

IF  (.NOT.  3ETC011)  GO  TO  210 

Ltine  invariant  calculations  of  ji] 
210  CONTINUE 

[time  varyinf;  calculations  of 
RETURN 
301  CONTINUE 

Ltime  varying  calculations  of  Hj 
RETURN 
401  CONTINUE 

IF  (.NOT.  SETCOli)  GO  TO  410 

[time  invariant  calculations  of  R] 
410  CONTINUE 

[time  varying  calculations  of  Hj 
RETURN 
END 


MSA 3.  Subroutine  I'ZkS  calculates  and  provides  to  ilCAP  the  noise- 
corrupted  system  measurements  using  either  Equation  (5)  or  (l2).  The 
subroutine  statement  contains  the  following  arguments:  X,  DUN,  Z,  N2, 
M,  T.  The  vector  X contains  the  ciurrent  values  of  the  system  states 
(an  NZ-dimensioned  vector).  DUM  is  an  N 2- dimensioned  vector  which 
can  be  used  by  the  subroutine  for  temporary  calculations.  The  calcu- 
lated measurements  are  returned  to  iIGAP  in  Z (an  N-dimensioned  vector). 
R2  is  the  dimension  of  the  system,  i'i  is  the  number  of  measurements 
of  the  system.  The  current  value  of  the  simulation  time  is  contained 
in  T.  Subroutine  NOISE  is  called  to  corrupt  the  system  measurements 
with  noise.  Table  VI  lists  the  required  Fortran  statements  for  sub- 
routine nsAS. 
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Table  VI 


i 

:| 


I 


Required  statements  for 

SUBROUTINE  lliiAS  (X,  DU.'I,  Z,  N2,  ll,  T) 

DlilEloION  X(l),  DUi:(l),  Z(l) 

[^measurement  calculations"] 

RETURN 

EI<D 

UPDATXF . Subroutine  UPDATXF  calculates  and  provides  to  JICAP  the 
updated  filter  estimate,  XF+,  The  subroutine  statement  contains  the 
following  arguments;  X,  FK,  Z,  Nl,  i!.  The  Kl-dimensioned  vector  X 
contains  the  current  filter  state  value  (XF-).  The  array  FK  (of 

t ] 

dimension  Klx;l)  contains  the  Kalman  gain  matrix  values.  The  noise  . 

system  measurements  are  contained  in  the  ll-dimensioned  vector  Z.  K1 
is  the  dimension  of  the  filter  and  li  the  number  of  measurements  of  the 
system.  The  updated  filter  estimate  (XF+)  is  returned  to  MCAP  in  X. 

UPDATXF  is  not  a constant  subroutine  because,  in  general,  nonlinear  ] 

calculations  may  be  necessary  Lsee  Equation  (20)3.  However,  a constant  I 

subroutine  IN-IXVE  is  available  to  perform  calculations  of  a matrix  | 

times  a vector,  2.  “ A.  arguments  for  KMXVE  are:  A,  3,  G,  N,  and  j 

M.  Where  A is  an  NxII  matrix,  B is  an  H-dimensional  vector,  and  G is  a j 

vector  of  dimension  JI  for  result.  Table  VII  lists  the  required  | 

Fortran  statements  for  subroutine  Ul^DATXF . ! 

I 

i 

\ 

1 

S 


Table  VII 

Required  Statements  for  UPDATXF 


SUBROUTINE  UPDATXF  (X,  FK,  Z,  Nl,  fl) 
DllkI.SICN  X(l),  FK(1),  Z(l) 

[^update  calculations^ 

RETURN 

END 
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GORRHCT.  Subroutine  GORRSGT  applies  impulsive  feedback  by  up- 
dating (selected)  system  states,  and  after  feedback,  resetting  the 
corresponding  filter  states  to  zero.  The  subroutine  statement  contains 
the  following  arguments:  X3,  XF,  112,  and  Nl,  The  K2-dimenoioned  vector 
XS  contains  the  current  values  of  the  system  states,  XF  is  an  Nl- 
dimensional  vector  that  contains  the  current  values  of  the  filter 
states.  N2  is  the  dimension  of  the  system.  Nl  is  the  dimension  of  the 
filter.  The  impulsive  feedback  is  performed  by  calculating  Equations 
(51)  and  (53)  and  returning  the  result  in  the  X3  and  XF  arrays. 

Table  VIII  lists  the  required  Fortran  statements  for  subroutine  GCRRilQT. 


U3RIH.  Subroutine  U3RIN  is  used  to  read  in  the  problem  dependent 
user-supplied  data  for  each  llonte  Carlo  analysis.  USRIN  may  be 
employ ed  in  .•>  variety  of  ways.  U3RIN  may  be  used  to  input  parameters 
which  change  between  Monte  Carlo  simulations  but  which  remain  unchanged 
within  a single  simulation.  Or  the  user  may  employ  USRIN  to  create 
special  initial  conditions.  • Whatever  the  usage,  U3RIK  is  a convenient 
method  which  permits  a user  to  customize  IlGAP  input  data  to  the 
specific  problem.  It  is  important  to  note  that  USRIN  is  called  only 
once  (during  function  Wl).  Table  IX  lists  the  required  Fortran  state- 
ments for  subroutine  USRIN . 
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Table  VIII 

Required  Statements  for  CORRECT 


SUBROUTINE  CORRECT  (X3,  XF,  N2,  Nl) 

(^impulsive  feedback  calculations^ 
RETUiu,' 

END 


Table  IX 

Required  3tatoments  for  U3UI1\ 


jU'JROUTlIlZ  U3RIII 

[user  specified  read^ 
RiTUKH 
'J;D 


3 

i 

i 

i 

I 


•! 

1 

i 

I 

i; 
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Section  V 
A Sample  Problem 

Introduction 

In  order  to  Illustrate  how  the  generalized  Honte  Carlo  Analysis 
Program  (flCAP)  can  be  used  to  analyze  and  evaluate  the  performance 
of  a suboptimal  Kalman  filter,  a sample  problem  is  illustrated  in 
this  section.  The  example  used  is  the  same  problem  investigated  in 
the  GGAP  users  guide  (Ref  2). 

The  problem  to  be  studied  is  a linear  system  from  which  multiple, 
noise-corrupted  measurem.ents  are  available.  The  truth  model  is 
presented  and  then  a suboptimal  Kalman  filter  is  designed  to  estimate 
certain  physical  (error)  states  of  this  system  precisely.  These 
models  are  then  programmed  into  the  appropriate  KCAP  subroutines  and 
10  runs  performed  for  a one  hour  performance  analysis. 

Problem  Formulation 

The  sample  problem  concerns  the  analysis  of  a single-ax.is  inertial 
navigation  system  shown  in  Figure  5.  In  the  figure,  R represents  the 
radius  of  the  earth  and  g is  the  acceleration  due  to  gravity.  The 
performance  of  the  system  is  analyzed  for  a one-hour  time  period, 
during  which  time,  position,  and  velocity  measurements  occur  at  30 
second  intervals.  The  integration  interval  is  30  seconds.  Printed 
output  is  required  every  40  integration  cycles  and  plotted  output  is 
required  every  integration  cycle. 

The  mathematical  model  of  the  inertial  navigation  system  errors 
is  formulated  using  state  variables.  That  is,  the  linear  time- 


! 1 


I 


j 
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varying  system  Is  described  by  the  vector  differential  equation* 

x(t)  - F(t)  x(t)  + G(t)  w(t)  (54) 

and  the  discrete  measurements  are  represented  byi 

v(t^)  (55) 

Considering  the  case  where  x(t),  w(t),  and  v(t)  are  Gaussian  random 
variables,  their  probability  density  functions  are  completely 
defined  if  their  first  and  second  moments  are  known.  The  following 
definitions  apply* 

x(t)  is  an  n-vector  denoting  the  system  state  with  zero 
meaji  and  covariance  P(t): 

ECi(t)  £(t)’']  - P(t) 

F(t)  is  the  nxn  plant  matrix 

w(t)  is  an  s vector  of  white  noise  inputs  which  corrupt 

the  states,  having  zero  mean  and  covariance  9t(t)6(T)* 

Stw(t)  w(t+T)'^]  = i(t)  6(t) 

G(t)  is  the  nxs  noise  matrix 

H(tj^)  is  the  mxn  measurement  matrix 

v(tj^)  is  an  ra  vector  of  white  noise  with  zero  mean  and 
covariance 

S[v(tj^)  v(tp  ] = R(t^)  6^^ 

Truth  Model.  The  single-axis  inertial  navigation  system  is 
described  by  the  following  differential  equations* 


(56) 

5^  B _g^  + b 

(57) 
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+ € 


(58) 


where  the  gyro  drift,  €,  is  modeled  as  the  sun  of  an  exponentially 
time  correlated  variable,  €^;  and  random  constant, 


. "1 


+ Wi 


(59) 


(60) 


where  the  accelerometer  bias  is  modeled  as  the  sum  of  two  distinct 


exponentially  time  correlated  variables,  b^  and  b2J 
b. 


b,  “ - -i—  + Wo 


(61) 


- + Wo 

'T'b2  3 


(62) 


and  where 

“ 3600  (seconds) 

“ 300  (seconds) 

'^b2  “ 3600  (seconds) 

Position  and  velocity  measurements  are  provided  to  the  filter. 
However,  the  system  measurements  are  corrupted  with  exponentially  time 
correlated  noises  AP  and  AV  in  addition  to  the  white  noises  vector  v, 
AP  and  AV  aure  modeled  asj 


AP 


AP 


(63) 


AV 


AV  . 

+ Wc 

T^V  5 


(64) 
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where t 


where I 


^ rad/sec 

“ 6.44  X 10  ^ ft/sec^ 

-2  2 

= 3*22  X 10  ft/sec 

0^  “ 3.0  X lo"*”^  ft 

-1  , 

O^y  “ 5.0  X 10  ft/sec 

a^p  “1.0  X lO"^^  ft 

O^y  “ 5.0  X 10  ^ ft/sec 

To  complete  the  specification  of  the  optimal  filter  model,  Table  X 
lists  the  initial  diagonal  values  of  the  apriori  covariance  matrix, 
with  the  off  diagonal  values  zero  initially. 

Table  X 


Covariance  Matrix  P(t  ) Elements 

' o' 


State  Variable 

J’iifV 

6P 

3.6  X lO’*’’^ 

6V 

+2 

1.0  X 10  ^ 

3.14  X 10"^ 

-I  *; 

2.35  X 10 

\ 

4.15  X lO"-^ 

^2 

-14 

5.88  X 10 

^2 

1.038X  10"^5 

AP 

9.0  X lO'^'^ 

AV 

2.5  X 10"^ 

1 
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Filter  Model . The  underlying  concept  for  the  development  of  the 
filter  model  is  to  delete  some  of  the  less  significant  or  non-dominant 
system  states.  Typically,  when  states  are  deleted  from  the  truth 
model,  the  noise  strengths  driving  the  model  are  increased  to  compen- 
sate for  the  deleted  states.  The  example  used  in  the  GCAP  users* 
guide  simply  deleted  states  bg,  AP,  and  AV  from  the  truth  model, 
without  any  increase  in  noise  strengths. 

The  deletion  of  these  states  yields  the  filter  state  vector 


- L6P  6V  i b^] 


and  measxirement  vector  z: 

— f m m"^ 

The  following  matrices  describe  the  filter  model t 


0 0 -g 

0 i 0 


0 0 0 - 


0 0 0 


0 0 0 1 0 
0 0 0 0 1 
0 0 0 0 0 
0 0 0 0 0 


0 0 0 0 0 


0 0 0 0 
10  0 0 


R 


f 


(75) 

(76) 


with* 


Sf  ~ £f  f 


0 0 0 
0 0 0 
0 0 0 

0 0 0 

0 0 0 


0 

0 


2a. 


2a 


0 

0 

0 

0 

2 


(77) 


With  these  a's  and  t's  having  the  same  values  as  for  the  truth  model. 
The  filter  initial  covariance  P(t^)  values  are  the  first  five  elements 
noted  in  Table  X, 

MCAP  Program  Im-plementat!l  on 

In  order  to  adapt  IICAP  to  this  example  problem,  the  user  must 
furnish  subroutines  USRIN,  XFDOT.  XSrOT.  TRAJ,  FLTiLiT.  NEAS.  UPDATXF. 
and  CORRECT,  Subroutine  USRIN  reads,  stores,  and  prints  the  constant 
parameters  betKeen  computer  runs,  XFDOT  calculates  the  derivatives 
of  the  filter  states,  XSDOT  calculates  the  derivatives  of  the  system 
states  and  adds  the  zero  mean  Gaussian  noise,  to  the  system  state 
values  at  the  end  of  the  integration  interval,  TRAJ  provides  the 
constants  or  time  varying  parameters  required  by  FLTUiAT , FLTiLAT  calcu' 
lates  the  nonzero  elements  of  the  filter  matrices  H,  and  R, 


263 
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HEAS  calculates  the  system  noise-corrupted  measurements,  UPDATXF 
calculates  the  updated  filter  estimate  XF+.  CORRECT  applies  impulsive 
feedback  by  updating  the  system  states  and,  after  feedback,  resets  the  ' 

filter  states,  i 

USRIN,  Subroutine  USRIN  is  usually  coded  first  in  order  to  specify 
the  labeled  COfMON  variables  used  in  the  other  user  supplied  subroutines. 

These  variables  are  those  parameters  which  may  vary  between  computer 
runs  but  remain  unchanged  within  a single  computer  run  (Monte  Carlo 
simulation) . Variables  typically  in  this  category  are  the  filter 
tuning  parameters*  CJaP*  this  problem.  Since 

the  user  may  also  want  to  analyze  the  sensitivity  of  the  filter  to 

changes  in  the  system  noise  strength  or  noise  correlation  time,  these  I 

parameters  are  also  read  in  by  USRIN . 

Figure  6 shows  the  USRIN  Fortran  coding  used  for  this,  example. 

These  parameters  are  input  and  printed  by  means  of  the  NAMELIST 
statement.  The  values  of  the  parameters  are  stored  in  the  labeled 

common  blocks:  QFNOIS,  RFNOIS,  TCFIL,  QSNOIS,  RSNCIS,  and  TGSYS  for  ' 

use  by  the  other  user  subroutines.  Table  XI  shows  the  relationship 

between  the  Fortran  variable  names  and  the  parameters  of  the  filter  - 

and  truth  model. 

XFDOT.  Figure  7 shows  the  XFDOT  Fortran  coding  used  for  this 
example.  Note  that  since  this  problem  models  error  states  and  full 
state  feedback,  the  values  for  the  filter  X and  XUOT  are  zero  for  all 

t.  Two  subroutine  constants  (RE  and  G)  are  used  in  the  first  section  i 

of  the  subroutine.  In  the  second  section,  note  that  the  derivative  | 
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SIIUK’OIIIINE  IJSFaN 


1H1S  IvOUMNE:  Fi'frtUS  rtNli  PRINTS  THE  USER  SUPPL lED  CONSTf^NTS . 

C0MMf.)N/0FNC)IS/0F(2)/RFN0IS/RP(2)/TCFIL/TAUF(2) 
C(JM(-triH/H!.!NOIS/C)S  ( S ) , QSO/RSNCIIS/R5  < 2 ) / TCSYS/TAUS  ( 5 ) 

NAMF:.L  1 ST/ 1 N/TAUF  f RF , RF 
NAMEl.  I S T/ 1 NN/T  AUS  r QS » CJSO , RS 


RFAIKUrlN) 

WRITE(6»IN) 

REATKSrINN) 

WRITE(6rlNN) 

RETURN 

END 

Figure  6.  Fortran  Coding  fcr  USRIN 


. SLItiRODTINE  XFUOT  ( X . XUOT  r M » IN TT  , NOI S » IJT » T IMF ) 
jK******#********#****)!'*******************:************#*******  ********* 

TFUS  R'OIITINE  CALFUI  ATES  THE  FILTER  DIFFERENTIAL  EOHATICINS 
XFDOT  ==  F(Xr  f l) 

**  NOTE  — SINCE  THIS  EXAMPLE  PROBLEM  MODELS  ERROR  STATES 

THEN  XF(T)  = 0.  FOR  ALL  TIME  IT  WOULD  BE  SUFFICIENT 
TO  USE 

RETURN 

END 


HOWEUFR  TO  SHOW  THE  METHOD  OF  CODING  THIS  ROUTINE  THE  FOLLOWING 
IS  PRESENTED 
D I MENS  I.  ON  X ( 1 ) T XDOT  ( 1 ) 

COMMON/TCEIL/TAHF (2) 

LOGICAL  INTTrNOIS 
IF  (.NOT.  INIT)  GO  TO  200 
C INITl  AL:r/ATION  OF  CONSTANTS 

RE  = 2.09E+7 
G = :t:>.2 
RETURN 

C XDOT  CALCULATIONS 

200  CONITNUE 

IE  (NOIS)  GO  TO  300 
XDOT(I)  = X(2)/RE 
XD0r<2)  = -X(3)*G  + X<S> 

XDOTCI)  = X(2)/RE  + X(4) 

. XDOT  (' 'll  = -X<4)/TAUF(1) 

XDOT(ri)  = -X<5)/TAUF(2) 

RETURN 

300  CONIINUE 
RETURN 
END 

Figure  7.  Fortran  Coding  for  XFDCT 
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Table  XI 


! 

I 

i 

I 


Correspondence  of  P'ortran  Variables  and  Model  Paxeuneters 


Fortran 

Model 

Filter 

System 

QF(1) 

QS  (1) 

^1 

QF(2) 

Q3(2) 

QS(3) 

^^2 

QS(4) 

^AP 

QS(5) 

^AV 

RF(l) 

RS(1) 

RF(2) 

R3(2) 

°'6V 

TAUF(l) 

TAU3(l) 

B 

TAUF(2) 

TAU3(2) 

TAU3(3) 

TAUS(4) 

^AP 

TAUS(5) 

^AV 

calculations  [F^(t)  Xf('t)D  are  directly  coded  from  the  filter  differ- 
ential equations.  The  values  for  the  time  correlated  vauriables  are 
obtained  from  the  'user  specified  labeled  COMMON,  Table  XII  shows  the 
correspondence  between  the  Fortran  variables  and  the  model  states, 
XSDOT,  Figure  8 shows  the  X3D0T  Fortran  coding  used  in  this 
example.  Subroutine  NOISE  is  used  to  obtain  a different  initial 
value  for  the  random  constant  Cg,  system  state  6,  for  each  run.  The 
second  section  performs  the  derivative  calculations  [Pg(t)  Xg(t)]]  for 
the  system.  Note  that  the  derivative  calculations  are  directly  coded 
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r.miRinir INI-  XEiimr  (XfXuuTrNr ini TrNOJS.nriTiME) 


****$$***  *M**:****1f.*t***************'******************  if  ************ 

THIS  ROUTINE  Crtl.CULAIES  THE.  SYSTEM  HI FE ERENT I Al,  EOUATIONS 
XSnOT  = F<XSfT)  t G+W 
I'lTMENSKIN  X(  1 ) rXrifJTI  1 ) 

COMMON/  I CSYS/I  AUS  ( Ei ) 

COMHUN/OSNO  T S/QS  < 5 ) . OSO 
lOGICAL  JNIIrNOIS 
IF  (.NtiT.  INIT)  GO  TO  200 
C INt  riAI..I7ATI0M  OF  CONSTANTS 

RE  = :!.09E+7 
G - S2.2 
XMEAN  = 0. 

CAM..  NOISE  (OBOrXMEANrXNOISE) 

X(A)  XNOISE 
RETURN 

c xnnr  cai.culations 

200  CONTINUE 

IF  (NOIS)  GO  TO  300 

XOOT(l)  = X<2)/RE 

XriOT(2)  = .-G*X(3)  + X(5)  + X(7) 

XIim<3)  = X(2)/RE  + X<.4)  + X<6)' 

XJinT(4)  = -X(4)/TAU!i(  1 ) 

XnOT(E^)  = -X<5)/TAUS(2) 

XMOKA)  =0. 

Xi:iOT(7)  = -X(7)/TAUS<3> 

XriOT(O)  = -X<8)/TAUS(4) 

XnOKV)  = -X(9)/TAUS<E-;> 

RETURN 

300  CONTINUE 

C XnOI  + NOISE  CALCULATIONS 

J=3 

MO  310  1=1 f5 

RMS  = OS(I)  * S0RT(2./TAUS< I ) > 

CALL  NOISE  (RMSrXMEAN. XNOISE) 

IF  (I  .EO.  3)  J=J+1 
J = J + 1 

X(J)  = X<J>  + XNOISE  * SORTdiT) 

310  CONTINUE 
RETURN 
ENM 

Figure  8,  Fortran  Coding  for  XSDOT 


SUMROUTINE  TRAJ(SETCON) 

THIS  ROUTINE  GENERATES  THE  FLIGHT  PROFILE. 

LOGICAL  SETCON 
COMMON/TRJCOM/RE»G 
IF < .N01 .SETCON)  GO  TO  100 

THE  FOLLOWING  CALCULATIONS  ARE  TIME  INVARIANT. 

RE=2.0VE+7 
0=32.2 
100  CONTINUE 
C 

C THE  FOLLOWING  CALCULATIONS  VARY  WITH  TIME. 

C 


RETURN 

ENM 


Figure  9,  Fortran  Coding  for  Traj 
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Table  XII 


ij 


[i 


Correspondence  of  Fortran  Vaxiables  and  Model  States 


Fortran 

Model 

State 

Filter 

System 

x(i) 

x(i) 

6P 

X(2) 

X(2) 

6V 

X(3) 

x(3) 

X(4)  • 

X(4) 

^1 

X(5) 

x(5) 

x(6) 

^2 

X(7) 

^2 

X(8) 

AP 

X(9) 

AV 

from  Equations  (56)  through  (62).  At  the  end  of  the  integration 
interval,  section  three  calculates  and  adds  it  to  the  system  state 
values.  In  this  example  the  only  noise  terms  are  those  added  to 
states  b^^,  b^,  AP,  and  AV,  which  correspond  to  state  numbers  4,  5» 
7,  8,  and  9. 

TRAJ.  Figure  9 shows  the  TRAJ  Fortran  coding  used  in  this  example. 
Due  to  the  simplicity  of  the  problem,  the  example  shown  is  trivial  but 
it  illustrates  how  TPAJ  might  be  used.  It  should  be  noted  that  the 
TRAJ  shoTO  here  is  identical  to  subroutine  TRAJ  used  in  the  GCAP  users 
guide  example. 
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FLTllAT.  Figure  10  shows  the  FLTilAT  Fortran  coding  used  in  this 
example.  Subroutine  FLTilAT  computes  the  nonzero  values  of  the  state 
space  matrices  for  the  filter  and  stores  these  values  in  the  variably 
dimensioned  array  A.  A separate  block  of  code  is  used  for  each 
matrix.  Within  each  block  of  coding,  the  nonzero  elements  of  that 
matrix  are  calculated  using  parameters  and  constants  transmitted  by 
labeled  COMMON  and  equated  to  the  appropriate  elements  of  the  A array. 

The  nonzero  elements  are  determined  by  examining  the  filter 
dynamics  matrix,  Equation  (73) • Row  1 column  2 is  the  first  nonzero 
element,  row  2 coluitm  3 is  the  second  nonzero  element,  and  so  on. 

Table  XIII  shows  the  correspondence  between  the  filter  matrices  non- 
zero elements  and  the  A array.  The  nonzero  values  for  matrices  F^, 

Hj.,  and  are  coded  using  Equations  (73) i (77) » (75) f anti  (76). 

It  should  be  noted  that  the  FLTllAT  shown  here  is  identical  to  subroutine 
FLTMAT  used  in  the  GCAP  users  guide  example, 

HEAS.  Figure  11  shows  the  HSAS  Fortran  coding  used  in  this  example. 
This  subroutine  calculates  and  provides  to  MCAP  the  noise-corrupted 
system  measurements.  The  system  measurements  Qig(t^)  Xs(ii)3  cal- 
culated and  the  results  of  subroutine  NOISE  added  to  obtain  the 
realized  value  of  the  measurement. 

UPDATXF.  Figure  12  shows  the  UPDATXF  Fortran  coding  used  in  this 
example.  This  subroutine  calculates  the  filter  estimate  XF+  by  solving 
the  filter  update  equation,  Equation  (8).  Subroutine  MMXVE  is  used 
to  multiply  the  filter  Kalman  gain  matrix  FK  times  the  measurement 
residual  f - li‘'(ii)2.f(ii)l  stored  in  variable  Z,  The  updated  filter 
estimate  XF+  is  returned  to  IlGAP  in  the  array  X. 
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SUliRCJI I r I HF.  FLTHrt  r ( Sm  CON  f MATfitl , A r IND ) 

C 

C THT.S  RIIUTINF  CAl.i:i)l  ArfiS)  THF  NllN-ZrRO  RLFMENTTk  OF  THE 

<-■ FltrER  MATRICES  FFr  OF.  HF  ANIi  RF.  STORIHli  THF:M  IN  THF  ARRAT  A 

C 

C— T-  IF  SFICON  .EO.  FALSE.  ONLY  TIME  VARIANT  ELEMENTS  ARE  OFNERATED 
C 

lORICAL  SCI  CON 
C 

COMMON/OFMO I S/OF  ( ? ) /RFNOI S/RF  ( .? ) /TCFIL/TAUF  < 2 ) 

COMriON/T  I MER/T I ME  . T I MUH . T I M INT 
COMMDN/IRJCOM/RE.G 
C 

niMFHSION  A(l).INIi(l) 

C 

C IIETERMINE  WHICH  MATRIX  IS  TO  BE  GENEr<ATED. 

C 

GO  10(101. 102. 103. 104). MATNO 

101  CONTINUE 

IF ( .NOT.SETCON)  GO  TO  1001 
C 

C CALCULATE  INVARIANT  ELEMENTS  OF  FF. 

C 

A(1)=1./RE 
A<2)=-G 
A<3)=1 . 

A(4)=1./RE 
A(5)-=l . 

A((i)--1./TAUF(1) 

.A(7)---l  ./TAUF<2> 

1001  CONTINUE 
C 

CALCULATE  TIME  VARIANT  ELEMENTS  OF  FF. 

RETURN 

102  CONTINUE 

IF( .NOT.SETCON)  GO  TO  1002 
C 

C CALCULAIE  INVARIANT  ELEMENTS  OF  OF. 

C 

• no  1022  11=1.2 

A ( 1 1 ) =2 . 0)(i0F  < 1 1 ) *0F  < 1 1 ) /TAUF  (II) 

1022  CUNIINUE 

1002  CONTINUE 
C 

c — --  CALCULATE  TIME  VARIANT  ELEMENTS  OF  OF. 

C 

RETURN 

103  CONTINUE 

IF( .NOT.SETCON)  GO  TO  1003 
C 

C CALCULATE  INVARIANT  ELEMENTS  OF  HF . 

C 

A ( 1 ) = 1 . 

A(2)=l . 

1003  CONTINUE 
C 

n CALCULATE  TIME  VARIANT  ELEMENTS  OF  HE. 

C 

RETURN 

104  CONTINUE 

IF ( .NOT.SETCON)  GO  TO  1004 
C 

C CALCULATE  INVARIANT  ELEMENTS  OF  RF . 

C 

A(1)=F(F(1)**2 
A(2)---RF  <2)»*2 

1004  CONTINUE 
C 

C •CALCULATE  TIME  VARIANT  ELEMENTS  OF  RF . 

C 

RETURN 

END 

Figure  10,  Fortran  Coding:  for  FLTilAT 

2?0 
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SUlilUKJI  INK  MKrt!)  (X»XIltn  FZrNrMf  D 
C 

C*****.**¥**t*******Wm****^******.*****iifilnli:******K*t***t)t.*t*.tl*ti*.t****$***: 

C 

C IHIS  liUUilNE  CALCDI.ATES  THE  N(:)ISE-CORRIJPrt:i:i  MEASlIREMEiJ  m USED 
C PY  THE  FILTER 

L'  Z = H(XSfT)  + y 

c 

COMHON/RSNOIli/Rfi  O ) 
niHENSIUN  X<1)»Z(1) 
xmt'an  = 0.  ■ 

CAI.I.  NOISE  (RS(  1 ) fXMCANrXNOISE) 

Z<1)  - X(l)  + X(8)  + XNOISE 
CALL  NOISE  <RS<2) fXMEANf XNOISE) 

Z(2)  = X(2)  + X(9>  + XNOISE 

RETURN 

END 

Figure  11,  -Fortran  Coding  for  MiiAS 
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SUPROIJTINE  URDATXF  (X,FK»ZrN»M) 

C 

C 

C THIS  ROUTINE  CALCULATES  XFT  GIUEN  THE  FILTER  GAIN  MATRIX  FK 

C AND  THE  SYSTEM  MEASUREMENTS f Z 

C 

C XFI  = XF-  -I  K*<Z(SYSTEM)  - HCXF-fT)) 

C 

COMMON  WRK(l) 

DIMENSION  X(1>fFK(1)fZ(1) 

DO  110  I“1fM 
ZCI)  =--  Z(I)  - X(I) 

110  CONTINUE 

CALL  MMXOE  ( FK f Z f URK f N f M ) 

DO  120  I--1fM 

X(I)  = X(I)  + WRK<I) 

120  CONTINUE 
RETURN 
END 

Figure  12.  Fortran  Coding  for  UPDATXF 


SLiiROiniNE  CORRECT  < XS  f XF  f IXSf  IXF ) 

« * * * « )!  >i!  ^ * )f:  K !l  >K  I:  t >1  !K  iK  K « )|(  !|(  )<(  !L  * « 4:  IK  « !K  liuL' Id  >1!  ^.  )|<  :i!  >1  * * >1 D:  !K  !K  IT  IT  4:  H iK  !|(  !|( « IT  « !K )!' Hc  1):  * He 

■ THIS  ROUTINE  APPLIES  INSTANTANEOUS  FEEDDACK  CORRECTION  TO 
SYSTEM  STATES  AND  RESETS  THE  FILTER  ESTIMATES 


DIMENSION  XS<1)fXF<1) 

DO  100  1=1 F IXF 
XS(I)  = XS(I)  - XF<I) 

XF<I)  = 0. 

100  CONTINUE 
. RETURN 
END 

Figure  13.  Fortran  Coding  for  CORHFCT 
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Table  XIII 


Correspondence  of  Nonzero  Elements  and  the  A Array 


Nonzero  Element 

Array 

A 

Matrix 

Row 

Column 

1 

2 

1 

2 

3 

2 

2 

5 

3 

^f 

3 

2 

4 

3 

4 

5 

4 

6 

5 

7 

B 

1 

H 

mi 

2 

H 

~f 

1 

1 

1 

2 

2 

2 

R 

“f 

1 

1 

1 

2 

2 

2 

CORRECT.  Figure  I3  shows  the  CORRECT  Fortran  coding  used  in  this 
example.  All  of  the  filter  states  estimavcsd  are  used  by  the  inertial 
navigation  system  to  zero  out  the  system  error  states.  Consequently, 
CORRECT  calculates  the  Instantaneous  feedback  correction  of  the  system 
states  using  Equation  (5l).  After  the  correction  has  been  applied,  the 


filter  states  are  reset  to  zero. 


fj 


!j 


■j 


HCAP  Sample  Input  Data 

Table  I listed  and  described  the  input  parameters  required  for 
HCAP.  Figure  14  shows  the  input  card  data  deck  used  in  this  example. 

As  mentioned  in  the  problem  formulation,  the  performance  of  the  iner* 
tial  navigation  system  is  to  be  analyzed  from  the  initial  time  (T^  = 0.) 
for  a one  hour  time  period  (TF  = 3^00.  seconds).  The  measurements  occur 
at  30  second  intervals  (DTUP  ■=  30.)  and  the  Integration  time  is  30 
seconds  (DTINT  = 30.).  Printed  output  is  required  every  40  integration 
cycles  (IPRGT  = 40)  for  two  of  the  simulation  runs  (IRNPRCT  =2).  The 
Kalman  gain  matrix  K is  also  printed  (IPRK  **  l).  Plotted  output  is 
required  every  integration  cycle  (IPLCT  = l).  Ten  runs  are  to  be  per- 
formed in  the  performance  analysis  (IPA3S  = lO).  The  random  number 
generator  is  started  from  a "seed"  of  77  (ISKED  = 77). 

The  number  of  nonzero  terms  and  indices  of  those  terms  is  obtained 
froi,;  Table  XIII  and  shown  in  the  proper  format  in  Figure  14.  Since 
error  states  are  being  modeled  and  since  this  is  an  example,  the  ini- 
tial conditions  for  the  filter  CXF(O)]]  and  system  [XS(0)]]  are  set  to 
zero.  The  filter  covariance  diagonal  element  initial  condition  C^F(0)3 
is  obtained  from  Table  X. 

The  U.SRIN  data  is  shown  in  Figure  14  as  the  NAHELIST  $IN  and  $INN 
group  names.  The  values  shown  in  the  figure  are  obtained  using  Table 
XI  and  the  problem  formulation. 

HCAP  Sample  Output  Data 

The  printed  data  from  the  Monte  Carlo  simulation  of  10  rxms  is 
shown  in  Figures  15  through  23.  The  printed  data  consists  of  the 
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MONTE  CARLO  ANALYSIS  PROGRAM  (MCAP)  EXAMPLE  CASE 

♦ USRCTL  TO  = O.Or  TF  = 3600.  Of  DTUP  = 30.  Or  IiTINT  = 30.f 
N1  = 5r  N2  = 9f  M = 2r 

IPASS  = JOf  IRNPRCT  = 2r  ISEED  = 77r 
IPRCT  = 40.  IIH.CT  -■-  1-r  IPRK  = 1 $ 

7 **  N/.rF  ** 

12  23  25  32  34  44  55 

2 **  N/OF  ** 

4 4 5 5 

2 **  N7IIF  ** 

112  2 
2 #*  NZRF  ** 

112  2 i 

0.  **  XF<0)  **  } 

0. 

J 

0.  . I 

0. 

0. 

0.  **  XS(0)  ** 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

3.A0OE+O7’  **  PF(0)  ** 

1.000Et02 
3.140E-04 
2.350E-15 
4. 150E-05 

$IN  TAUFCl)  = 2*3600. f 
0F(1)  = 4.05E-8f  6.44E-3f 

Pr<l)  “ IOC. Or  0.5r  $ 

♦ INN  TAUSd)  = 2*3600..  300.0.  1800.0.  300.0. 

QS(1)  = 4.85E-8.  6.44E-3.  3.22E-2.  3.0E+2.  5.0E-1. 

RS<1>  = 100.0.  0.5.  * 


Figure  14.  Data  Input  Card  Deck 


tUS^C'V. 


CfN^pst.  •40NTE  C»i?LO  ANALYSIS  FRCGRAM  12/09/77  17.15.02.  RUN  NUMBER 

MONTr  CARLO  ANALYSIS  PPOGPAM  (MCAP)  EXAMPLE  CASE 
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following  Items t 


1.  The  problem  title,  date,  and  time  the  simulation  was  conducted. 

2.  The  echo  printing  of  $USRCTL. 

3.  The  echo  printing  of  the  nmber  of  nonzero  elements  for  each 
state  space  matrix  together  with  the  corresponding  indices  of 
these  elements. 

4.  The  minimum  required  dimension  of  the  unlabeled  COillON  block. 

5.  The  echo  printing  of  U3RIN  data  $IN  amd  $INH. 

6.  The  initial  values  of  the  filter  and  system  states  []XF(o)  and 
XS(0)]. 

7.  The  square  root  of  the  initial  diagonal  values  of  the  filter 
covariance  matrix 

8.  The  values  of  the  filter  and  system  states  before  and  after  the 
measurement  and  after  the  impulsive  correction, 

9.  The  values  of  the  measurements, 

10,  The  square  root  of  the  diagonal  elements  of  the  filter  covari- 
ance matrix  at  times  other  than  and  TF, 

11,  The  full  Kalman  gain  matrix,  K,  at  the  update  times. 

12,  The  final  lower  triangle  values  of  the  filter  covariance  matrix. 


written  onto  TAPii3  UCAP  and  then  used  hy  a separate  Calcomp 
plotting  program.  Figures  24  through  28  are  the  standard  deviation 
plots  and  Figures  29  through  33  are  the  mean  error  + s plots  for 
this  example.  It  should  be  noted  that  this  example  has  not  been 
tuned. 


FI&  27. 

LINERR  FILTER  EXRMPLE  (10  RUNS) 

287 

1 



1 

1 

. ..  J 

00*8 

oo*v 

00*0 

OO't'-  00*8-° 

A3G 

019  -/+ 

Nb3W  (3)X 

FIG  31 


LINERR  FILTER  EXRMPLE  (10  RUNS) 

291 


so. 00  160.00  240.00  320.00  400.00  480.00  660.00 

TIME  (SEC)  *10' 


o 

o 


o 


FIO  33.  LINERR  FILTER  EXPHPLE  (10  RUNS) 

293 


Appendix  B 


Linearization  of  the  Sub-Optimal  Filter  State  and  lleasurement 


Equations 


To  develop  the  filter  model,  the  first  six  states  of  the  system 
model  were  deleted  and  the  inertial  acceleration  of  the  satellite  with 
respect  to  the  tracker  in  the  tracker  frame,  (^g),  was  evaluated  as 
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